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Abstract 


This  document  presents  advanced  spectral  estimation  techniques  for  multiple  moving 
target  feature  extraction,  spectral  analysis  of  one-dimensional  (1-D)  gapped  data  sequences, 
and  synthetic  aperture  radar  (SAR)  imaging  with  angle  diversity  data  fusion. 

For  the  multiple  moving  target  feature  extraction  problem,  we  study  clutter  suppression 
and  feature  extraction  of  multiple  moving  targets  for  airborne  high  range  resolution  (HRR) 
phased  array  radar.  We  show  how  to  use  a  Vector  Auto- Regressive  (VAR)  filtering  tech¬ 
nique  to  suppress  the  correlated  ground  clutter  and  propose  a  relaxation-based  parameter 
estimation  algorithm  for  multiple  moving  target  feature  extraction.  For  the  spectral  analy¬ 
sis  of  gapped  data  sequences,  we  present  an  algorithm  for  nonparametric  complex  spectral 
estimation  of  gapped  data  via  an  adaptive  filtering  approach,  referred  to  as  the  GAPES 
(Gapped-data  Amplitude  and  Phase  Estimation)  algorithm,  which  iterates  between  step- 
s  of  estimating  the  adaptive  filter  and  the  corresponding  spectrum  via  APES  (Amplitude 
and  Phase  Estimation)  and  filling  in  the  gaps  via  a  least  squares  (LS)  fitting.  For  the  SAR 
imaging  with  angle  diversity  data  fusion,  we  study  both  nonparametric  and  quasi-parametric 
methods.  For  the  nonparametric  method,  we  propose  to  use  the  GAPES  algorithm  to  fill 
in  the  gaps  in  a  range-wise  mode  and  then  use  a  two-dimensional  (2-D)  APES  algorithm 
to  obtain  the  final  2-D  SAR  image.  For  the  quasi-parametric  method,  we  first  establish 
a  flexible  data  model  that  describes  each  target  scatterer  as  a  2-D  complex  sequence  with 
arbitrary  amplitude  and  constant  phase  in  range  and  cross-range.  Then  we  present  an  algo¬ 
rithm,  referred  to  as  QUALE  (QUasi-parametric  ALgorithm  for  target  feature  Extraction), 
for  the  SAR  target  feature  extraction  and  imaging  with  angle  diversity  data  fusion  based  on 
the  flexible  data  model.  Numerical  results  are  presented  to  illustrate  the  performances  of  all 
algorithms  proposed  in  the  document. 


1.  Introduction 


This  document  presents  advanced  spectral  estimation  techniques  for  multiple  moving 
target  feature  extraction,  spectral  analysis  of  one-dimensional  (1-D)  gapped  data  sequences, 
and  synthetic  aperture  radar  (SAR)  imaging  with  angle  diversity  data  fusion. 

For  the  multiple  moving  target  feature  extraction  problem,  we  study  clutter  suppression 
and  feature  extraction  of  multiple  moving  targets  for  airborne  high  range  resolution  (HRR) 
phased  array  radar.  We  show  how  to  use  a  Vector  Auto-Regressive  (VAR)  filtering  tech¬ 
nique  to  suppress  the  correlated  ground  clutter  and  propose  a  relaxation-based  parameter 
estimation  algorithm  for  multiple  moving  target  feature  extraction.  For  the  spectral  analy¬ 
sis  of  gapped  data  sequences,  we  present  an  algorithm  for  nonparametric  complex  spectral 
estimation  of  gapped  data  via  an  adaptive  filtering  approach,  referred  to  as  the  GAPES 
(Gapped-data  Amplitude  and  Phase  Estimation)  algorithm,  which  iterates  between  step- 
s  of  estimating  the  adaptive  filter  and  the  corresponding  spectrum  via  APES  (Amplitude 
and  Phase  Estimation)  and  filling  in  the  gaps  via  a  least  squares  (LS)  fitting.  For  the  SAR 
imaging  with  angle  diversity  data  fusion,  we  study  both  nonparametric  and  quasi-parametric 
methods.  For  the  nonparametric  method,  we  propose  to  use  the  GAPES  algorithm  to  fill 
in  the  gaps  in  a  range-wise  mode  and  then  use  a  two-dimensional  (2-D)  APES  algorithm 
to  obtain  the  final  2-D  SAR  image.  For  the  quasi-parametric  method,  we  first  establish 
a  flexible  data  model  that  describes  each  target  scatterer  as  a  2-D  complex  sequence  with 
arbitrary  amplitude  and  constant  phase  in  range  and  cross-range.  Then  we  present  an  algo¬ 
rithm,  referred  to  as  QUALE  (QUasi-parametric  ALgorithm  for  target  feature  Extraction), 
for  the  SAR  target  feature  extraction  and  imaging  with  angle  diversity  data  fusion  based  on 
the  flexible  data  model. 

This  document  contains  5  chapters.  In  Chapter  2,  we  study  the  multiple  moving  target 
feature  extraction  for  airborne  HRR  phased  array  radar.  Chapter  3  presents  a  technique  for 
nonparametric  spectral  analysis  of  gapped  data  sequences.  In  Chapters  4  and  5,  we  study 
the  SAR  imaging  with  angle  diversity  data  fusion  via  nonparametric  and  quasi-parametric 
methods,  respectively. 
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In  Chapter  2,  we  study  the  clutter  suppression  and  feature  extraction  of  multiple  moving 
targets  for  HRR  phased  array  radar.  To  avoid  range  migration  problems  that  occur  in  the 
HRR  radar  data,  we  divide  each  HRR  profile  into  non-overlapping  low  range  resolution 
(LRR)  segments,  so  that  each  LRR  segment  contains  a  sequence  of  HRR  range  bins.  No 
information  is  lost  due  to  the  division  and  hence  no  loss  of  resolution  occurs.  We  show 
how  to  use  the  VAR  filtering  technique  to  suppress  the  correlated  ground  clutter.  Then 
a  relaxation-based  parameter  estimation  algorithm  is  presented  for  multiple  moving  target 
feature  extraction.  The  problem  of  multiple  target  feature  extraction  is  reduced  to  the  feature 
extraction  of  a  single  target  in  a  relaxation-based  iteration  step.  In  each  iteration  step  and 
for  each  target,  the  target  phase  history  sequence  and  Direction-Of-Arrival  (DOA)  (or  the 
unknown  array  manifold)  are  estimated  from  some  spatial  signature  vectors  by  minimizing 
a  Weighted  Least  Squares  (WLS)  cost  function.  The  complex  amplitude  and  range  of  each 
target  scatterer  are  then  extracted  from  the  estimated  target  phase  history  sequence  by  using 
RELAX,  a  target  feature  extraction  algorithm  with  super  resolution  performance. 

In  Chapter  3,  We  present  the  GAPES  algorithm  for  nonparametric  complex  spectral 
estimation  of  gapped  data.  Unlike  the  fast  Fourier  transform  (EFT),  windowed  or  averaged 
EFT  spectra,  the  APES  spectrum  has  good  resolution  properties,  suffers  from  little  or  no 
leakage  effects,  and  has  good  statistical  stability.  The  excellent  performance  of  APES  in 
this  class  of  nonparametric  spectral  analysis  methods  is  one  of  the  reasons  why  we  choose 
to  extend  this  particular  approach  to  the  gapped-data  case.  The  incomplete  data  sequence 
may  contain  gaps  of  various  sizes.  The  GAPES  algorithm  iterates  the  following  two  steps; 
(1)  estimating  the  adaptive  FIR  filter  and  the  corresponding  complex  spectrum  via  APES, 
and  (2)  filling  in  the  gaps  via  a  LS  fitting  criterion.  The  initial  condition  for  the  iteration  is 
obtained  from  the  available  data  segments  via  APES. 

In  Chapter  4,  we  investigate  the  SAR  imaging  with  angle  diversity  data  fusion  via 
GAPES.  The  APES  and  GAPES  algorithms  are  first  introduced  for  the  spectral  estima¬ 
tions  of  complete  and  gapped  data,  respectively.  For  the  angle  diversity  data  fusion  for  SAR 
imaging,  the  radar  measurements  in  range  are  complete  whereas  the  gaps  are  caused  by  the 
intermittent  measurement  mode  and  thus  only  exist  in  the  azimuth  (cross-range)  dimen- 
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sion.  We  perform  the  1-D  windowed  FFTs  (WFFTs)  in  range,  use  the  GAPES  algorithm 
to  interpolate  the  gaps  in  the  aperture  for  each  range,  apply  the  1-D  inverse  FFTs  (IFFT- 
s)  and  de-window  in  range,  and  apply  the  two-dimensional  (2-D)  APES  algorithm  to  the 
interpolated  matrix  to  obtain  the  final  2-D  SAR  image. 

In  Chapter  5,  We  first  establish  a  flexible  data  model  that  describes  each  target  scatterer 
as  a  2-D  complex  sinusoid  with  arbitrary  amplitude  and  constant  phase  in  range  and  cross¬ 
range,  and  then  present  the  QUALE  algorithm  for  SAR  target  feature  extraction  and  imaging 
via  data  fusion  through  angle  diversity  based  on  the  established  flexible  data  model.  QUALE 
first  estimates  the  model  parameters  that  include,  for  each  scatterer,  a  2-D  arbitrary  real¬ 
valued  amplitude  sequence,  a  constant  phase,  and  scatterer  locations  in  range  and  cross- 
range.  QUALE  then  averages  the  estimated  2-D  real-valued  amplitude  sequence  over  range 
by  making  the  assumption  that  the  scatterer  radar  cross  section  (RCS)  is  approximately 
constant.  QUALE  next  models  the  so-obtained  1-D  sequence  with  a  simple  sine  function 
by  assuming  that  the  scatterer  is  approximately  a  dihedral  (a  trihedral  is  approximated  as 
a  very  short  dihedral) ,  and  estimates  the  relevant  sine  function  parameters  by  minimizing  a 
nonlinear  least  squares  (NLS)  fitting  function.  Finally,  the  approximate  2-D  SAR  image  is 
reconstructed  by  using  the  estimated  features. 

Each  of  the  aforementioned  chapters  is  self-contained  with  its  own  introduction,  for¬ 
mulation  of  the  problem  of  interest,  detailed  presentation  of  approaches,  conclusions,  and 
references.  We  acknowledge  Dr.  Rob  Williams  of  Air  Force  Research  Laboratory  for  giving 
us  the  interesting  research  topic  of  SAR  imaging  with  angle  diversity  data  fusion. 
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2.  Multiple  Moving  Target  Feature  Extraction  for 

Airborne  HRR  Radar 


2.1  Introduction 

Clutter  and  jamming  suppression  is  critical  for  airborne  radar  signal  processing.  The 
ground  clutter  observed  by  an  airborne  radar  is  spread  over  two  dimensions  of  both  the  range 
and  spatial  angle  and  the  clutter  spectrum  also  covers  a  certain  Doppler  region  due  to  the 
platform  motion.  A  Vector  Auto-Regressive  (VAR)  filtering  technique  was  recently  proposed 
by  Swindlehurst  and  Stoica  [1]  to  suppress  the  clutter  adaptively.  It  is  more  robust  than 
the  non-adaptive  Displaced-Phase-Center-Antenna  (DPCA)  processing  since  the  latter  is 
sensitive  to  both  array  calibration  errors  and  system  mismatch.  The  VAR  filtering  technique 
whitens  the  correlated  clutter  only  temporally  and  can  be  computationally  simpler  than  the 
conventional  Space-Time  Adaptive  Processing  (STAP)  based  techniques,  which  may  require 
a  significant  amount  of  computations  due  to  the  need  of  using  a  bank  of  filters  and  the 
inversions  of  matrices  of  large  dimensions,  STAP  cannot  be  used  for  jamming  suppression 
if  the  secondary  data  bins  used  to  obtain  the  second-order  statistics  of  the  ground  clutter 
do  not  contain  jamming  interference.  Although  the  VAR  filtering  technique  can  be  easily 
used  for  spatial  whitening  as  well,  it  is  not  needed  since  the  VAR-filtered  interference  is 
assumed  to  be  spatially  colored  with  an  unknown  and  arbitrary  covariance  matrix,  which 
automatically  achieves  jamming  suppression  when  the  VAR  filter  output  is  used  with  the 
Maximum  Likelihood  (ML)  methods  presented  by  Swindlehurst  and  Stoica  [1]  to  estimate 
the  target  parameters  for  Low  Range  Resolution  (LRR)  wide-area  surveillance  radar. 

Future  airborne  radars  will  be  required  to  provide  increasingly  High  Range  Resolution 
(HRR)  features  of  ground  targets,  which  makes  the  signal  processing  needed  by  airborne 
HRR  phased  array  radars  even  more  important.  Compared  to  a  conventional  airborne  LRR 
radar  [2],  an  airborne  HRR  radar  can  not  only  enhance  the  radar’s  capability  of  detecting, 
locating,  and  tracking  moving  targets,  but  also  can  provide  more  features  for  applications 
including  Automatic  Target  Recognition  (ATR)  [3, 4].  To  avoid  the  range  migration  problems 
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that  occur  in  HRR  radar  data,  we  first  divide  the  HRR  range  profiles  into  LRR  segments. 
Since  each  LRR  segment  contains  a  sequence  of  HRR  range  bins,  no  information  is  lost 
due  to  the  division  and  hence  no  loss  of  resolution  occurs.  We  extend  the  VAR  filtering 
technique  and  the  unstructured  ML  method  in  [1]  for  clutter  suppression  and  single  target 
feature  extraction  for  airborne  HRR  phased  array  radar.  The  target  phase  history  sequence 
and  Direction-Of- Arrival  (DOA)  (or  the  unknown  array  manifold)  are  estimated  from  some 
spatial  signature  vectors  by  minimizing  a  Weighted  Least  Squares  (WLS)  cost  function,  and 
the  complex  amplitude  and  range  of  each  target  scatterer  are  extracted  from  the  estimated 
target  phase  history  sequence  by  using  RELAX  [5]. 

We  also  extend  the  single  target  approaches  to  the  case  of  multiple  targets.  The  multi¬ 
ple  moving  target  scenario  occurs  frequently  in  radar  applications.  Yet  to  the  best  of  our 
knowledge,  little  research  on  the  topic  has  been  reported  in  the  literature.  We  present  a 
relaxation-based  algorithm  for  multiple  moving  target  feature  extraction.  Each  of  the  tar¬ 
gets  is  assumed  to  have  a  rigid-body  and  the  scatterers  of  the  same  target  have  the  same 
DOA.  The  relaxation-based  algorithm  is  used  to  minimize  a  nonlinear  least  squares  fitting 
function  by  letting  only  the  parameters  of  one  target  vary  at  a  time  while  fixing  the  para¬ 
meters  of  all  other  targets  at  their  most  recently  determined  values.  Thus  the  problem  of 
multiple  target  feature  extraction  is  reduced  to  the  feature  extraction  of  a  single  target  in  a 
relaxation-based  iteration  step.  We  use  numerical  examples  to  demonstrate  the  performance 
of  this  algorithm  for  clutter  suppression  and  multiple  moving  target  feature  extraction. 

The  remainder  of  this  paper  is  organized  as  follows.  In  Section  2,  we  establish  the  multiple 
moving  target  data  model  for  airborne  HRR  phased  array  radar,  which  is  followed  by  a  brief 
discussion  of  the  VAR  filtering  technique.  In  Section  3,  we  present  the  relaxation-based 
multiple  moving  target  feature  extraction  method.  Simulation  results  and  their  analysis  are 
presented  in  Section  4.  Finally,  we  give  the  conclusions  in  Section  5. 
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2.2  Data  Model  and  VAR  Filtering 


The  range  resolution  of  a  radar  is  determined  by  the  transmitted  signal  bandwidth.  To 
achieve  high  range  resolution,  a  radar  must  transmit  wide  band  pulses,  which  are  often  linear 
frequency  modulated  (LFM)  chirp  signals  [6].  The  range  resolution  of  a  LRR  radar  is  much 
larger  than  the  length  of  a  target  so  that  the  target  occupies  only  one  LRR  range  bin,  while 
the  range  resolution  of  a  HRR  radar  is  so  small  that  each  target  occupies  several  HRR  range 
bins. 

Consider  an  airborne  HRR  radar  having  a  one-dimensional  (1-D)  antenna  array  with  M 
elements  uniformly  spaced  along  the  flight  path  of  an  airborne  platform.  A  cluster  of  N  chirp 
pulses  is  transmitted  during  a  coherent  processing  interval  (CPI).  After  dechirping,  sampling, 
and  Fourier  transforming  the  signals  in  each  element  of  the  array,  we  obtain  the  HRR  range 
profiles.  Without  clutter  and  jamming  suppression,  these  clutter  and  jamming  dominated 
profiles  are  not  useful  for  any  applications.  To  avoid  range  migration  problems,  we  divide 
each  HRR  profile  into  non-overlapping  LRR  segments  so  that  each  LRR  segment  contains  L 
HRR  range  bins,  as  shown  in  Figure  2.1.  We  choose  L  to  be  much  larger  than  the  maximum 
number  of  range  bins  over  which  all  targets  can  possibly  expand  and  migrate  during  the 
CPI.  We  then  apply  inverse  Fourier  transform  to  each  segment.  For  the  segment  of  interest, 
where  targets  may  be  present,  the  inverse  Fourier  transform  yields  the  primary  data.  We 
assume  that  D  targets  are  present  in  the  primary  data  with  the  dth  target  consisting  of 
scatterers.  We  assume  that  the  scatterers  of  each  target  have  the  same  Doppler  frequency 
and  the  same  DOA,  but  different  complex  amplitudes  and  range  frequencies.  Then  the 
primary  data  model  can  be  written  as  (see  Appendix  A  for  the  model  derivation): 

D  (K^  \ 

Xi(n)  =  E  E  +  ei(n), 

\k=i  J 

Z  =  0,---,L-1,  n  =  0,-'-,iV-l,  (2.1) 

where  Xi(n)  is  the  array  output  vector  of  the  Zth  phase  history  sample  due  to  the  nth  pulse; 
a(0d)  is  the  array  manifold  and  is  a  function  of  the  dth  target  DOA  6d  relative  to  the  flight 
path;  ei{n)  is  the  interference  including  the  temporally  and  spatially  correlated  Gaussian 
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ground  clutter,  both  temporally  and  spatially  white  Gaussian  noise,  and  possibly  a  jammer 
that  is  temporally  white  but  a  point  source  in  space.  We  assume  that  the  clutter,  noise, 
and  jamming  in  diflferent  HRR  range  bins  are  independent  and  identically  distributed.  The 
complex  amplitude  Odk  and  the  frequency  fdk  are,  respectively,  proportional  to  the  RCS  and 
range  of  the  kth  scatterer  of  the  dth  target.  The  Vd  and  Ud  are,  respectively,  the  scaled  radial 
velocity  and  the  normalized  Doppler  frequency  of  the  dth  target.  Range  migration  can  occur 
due  to  the  radial  motion  between  the  radar  and  target  and  the  high  range  resolution  of  the 
HRR  radar.  For  notational  convenience,  let 

ujca  =  u}d{'i-  +  rl),  i  =  0,  d  =  l,  •••,!),  (2.2) 

where  r  =  is  a  known  constant  independent  of  the  target  motion  (see  Appendix  A) 

and  is  usually  very  small  (<C  0.01).  Then  (2.1)  can  be  written  as 

D 

XiW  =  +  l  =  0,---,L-l,  n  =  0,---,iV- 1,  (2.3) 

d=l 

where  we  drop  the  dependence  of  ad  on  9d  for  notational  brevity  and 
Note  that,  when  L  =  1,  the  model  in  (2.3)  reduces  to  the  data  model  for  the  LRR  case. 
For  L  >  1,  we  have  a  phase  history  sequence  for  each  LRR  segment  and  no  loss  of  range 
resolution  occurs  because  of  no  information  loss. 

The  secondary  data  are  obtained  from  segments  adjacent  to  the  segment  of  interest  in 
the  same  way  as  the  primary  data  are  obtained  from  the  segment  of  interest  (see  Figure  2.1. 
The  secondary  data  are  assumed  to  be  target  free  (see  [7]  for  the  latest  research  results  on 
the  selection  of  secondary  data)  and  are  modeled  as  a  VAR  random  process  [1].  The  VAR 
filter  has  the  form: 

+  (2.4) 

p=l 

where  P  is  the  VAR  filter  order,  z~^  denotes  the  unit  delay  operator,  and  I  is  the  identity 
matrix. 

Next,  we  outline  the  estimation  of  the  VAR  filter  using  the  target-free  secondary 

data.  Let  esi(n),  s  =  1,  •  •  • ,  S,  1  =  0,  •  •  • ,  L—1,  n  =  0,  •  •  • ,  —  1,  denote  the  secondary 
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data,  where  S  denotes  the  number  of  secondary  data  segments.  PVom 


L-l  s  N-\ 


Hi,---,Hp  =  arg  min 


with  II 'll  denoting  the  Euclidean  norm,  we  get 


6s/ (^)  +  X^Hpesj(n  p)||  , 

p=i 


(2.5) 


H  =  \  (2.6) 

where 

=  ['i'lo  •  •  •  ^sl  •  •  •  ’®^S(L-l)], 

E  =  [Elo  Ejj  •••  E5(£,_i)], 

'^si  =  bl’si{P)  ••• 

V’s/(’^)  =  -  [e]i(n  -  1)  •  •  •  ef,(n  -  P)]^ , 

and 

E,j  =  [e,i(P)  •••  e,,(Ar-l)],  Z  =  0,---,L-1,  s  =  l,---,P, 

with  (•)^  and  (•)^  denoting  the  transpose  and  conjugate  transpose,  respectively. 

Once  the  VAR  filter  coefficients  are  determined,  we  use  the  filter  to  suppress  the  clutter 
in  the  primary  data.  The  VAR  filter  output  for  the  primary  data  has  the  form 


y/(n)  =  il{z  ^)x,(n) 

=  ^’d/'Wd/a</e^“‘""  +  €/(n), 

<i=l 

Z  =  0,---,L-1,  n  =  P,-*-,N-l,  (2.7) 

where  'H{z~^)  has  the  same  form  as  'H{z~^)  in  (2.4)  except  that  the  {Hp}p_i  in  (2.4)  are 
replaced  by  {Hp}JLi, 

^dz  =  I  +  f:Hpe-™,  1  =  0,---,L-1,  d=l,. ••,£>,  (2.8) 

p=i 

and 


e,(n)  =  'k{z-^)ei{n),  Z  =  0,---,L-1,  n  =  P,  •  •  • ,  AT  -  1.  (2.9) 
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Let 


(^dl  —  /  =  0,  •  •  •  ,  jL  —  1,  d=l,---,D,  (2.10) 

be  referred  to  as  the  spatial  signature  vector  of  the  dth  target  for  the  /th  phase  history  sample. 
Then  (2.7)  can  be  rewritten  as 

yi{n)  =  E  +  e^n),  i  =  0,  •  •  • , L  -  1,  n  =  F,  •  •  • , iV  -  1.  (2.11) 

d=l 

Our  problem  of  interest  is  to  estimate  9d,  {o'dfc,  fdk]k=i}2=i  if  3-rray  manifold 
{a{9d)}d=i  is  known  or  {ud,  ad,  {odk,  fdk}k=i}d=i  if  ^-^ray  manifold  {a(0d)}f^i  is  un¬ 
known  from  the  VAR  filter  output  yi{n),  I  =  0,  -  •  ■  ,L-1,  n  =  P,---,  N-1,  by  minimizing 
Nonlinear  Least  Squares  (NLS)  criteria  using  a  relaxation-based  optimization  algorithm. 

2.3  Feature  Extraction  of  Multiple  Moving  Targets 

Our  feature  extraction  algorithm  consists  of  the  following  two  separate  steps: 

Step  I:  Estimate  the  target  space-time  parameters  9d,  if  the  array  mani¬ 
fold  {a(0d)}^i  is  known  or  {cOd,  a^,  if  the  array  manifold  is  unknown  from  the 

VAR  filter  output  {yi(n)},  /  =  0,  •  •  • , L  —  1,  n  =  P,  -  •  •  ,N  —  1, 

Step  II:  Estimate  the  target  range  parameters  {adk,  fdk}k=i^  from  the  esti¬ 
mate  of  {bdi}i'=Q  obtained  in  Step  I. 

2.3.1  Space-Time  Parameter  Estimation 

The  NLS  estimates  of  the  space-time  parameters  {oJd,  dd,  or 

{(jjd,  ad,  {bdi}i':^o}d=i  obtained  by  minimizing  the  following  NLS  criterion: 

L-l  D  2 

Ci  =  E  •  (2-12) 

1=0  d=l 

The  minimization  of  the  cost  function  Ci  in  (2.12)  is  a  highly  nonlinear  complicated  opti¬ 
mization  problem.  Here  we  present  an  alternating  or  relaxation-based  optimization  approach, 
which  is  a  conceptually  and  computationally  simple  method  for  multiple  moving  target  fea¬ 
ture  extraction.  The  relaxation-based  algorithm  is  used  to  minimize  Ci  by  letting  only  the 
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parameters  of  one  target  vary  at  a  time  while  fixing  the  parameters  of  all  other  targets  at 
their  most  recently  determined  values.  Therefore,  the  feature  extraction  of  multiple  moving 
targets  is  reduced  to  the  feature  extraction  of  a  single  moving  target  in  a  relaxation-based 
iteration  step.  We  first  consider  the  space-time  parameter  estimation  of  the  dth  target  and 
then  give  detailed  steps  of  our  approach  for  multiple  targets. 

Space-Time  Parameter  Estimation  of  the  dth  Target 
Let 

D 

yAn)=yiin)-  E  l  =  n  =  P,  •  •  • ,  AT  -  1,  (2.13) 

where  {a)j,  a^,  is  assumed  available.  Note  that  if  the  array  manifold  is 

known,  is  replaced  by  a(0i),  with  assumed  available.  Hence  yA''^) 

written  as 

ydi(n)  =  ad<e’‘^‘"”  +  edj(n),  l  =  0,---,L-l,  n  =  P,  •  •  •  ,iV  -  1,  (2.14) 

where  edj(n)  denotes  the  interference  due  to  clutter,  noise,  and  contributions  from  other 
targets.  We  assume  that  {€dj(n)}  is  a  zero-mean  temporally  white  Gaussian  random  process 
with  an  unknown  arbitrary  spatial  covariance  matrix  Qd.  Then  the  negative  log-likelihood 
function  of  yAA  (2-14),  is  proportional  to 

C2  =  ln|Qd|+Tr(Qj'Cd),  (2.15) 

where  [  •  |  and  Tr(-)  denote,  respectively,  the  determinant  and  the  trace  of  a  matrix  and 

Crf  =  E  X)  [y<«(^)  -  [yA'A  - 

i=0  Tl=:P 

=  (2.16) 
1=0 

with 

Yj,  =  yji(A'-l)l.  i  =  0,  --,L-l,  (2.17) 
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and 


l^di  =  ,  Z  =  0,  •  •  • ,  L  -  1.  (2.18) 

Minimizing  C2  in  (2.15)  with  respect  to  Qd  results  in  Qj  =  C^.  The  cost  function  in  (2.15) 
with  Qd  replaced  by  Qd  becomes: 


C,  = 


ln|Cd 
In 


L-l 

E 

NLto 


N  adi  — 


_  Ydi/3di\ 

j 


H 


+  Y. 


dl^dl  ~ 


Y.fl/3.„/3SYS 


N 


(2.19) 


The  minimization  of  C3  in  (2.19)  with  respect  to  oidi  gives 


Odi  =  -^YdiPdi 


Z  =  0,---,L-1. 


(2.20) 


The  estimate  of  the  Doppler  frequency  Ud  of  the  dth  target  is  obtained  by  minimizing  the 
concentrated  cost  function, 


074  =  In 


1 

m, 


L-l 


E 


Y^Y 


H 

dl  “ 


Y^P^I3«Y«- 

N 


(2.21) 


which  requires  a  1-D  search.  Note  that  the  spatial  signature  estimate  in  (2.20)  can  be 
interpreted  as  the  temporal  average  of  the  Doppler  shift  and  range  migration  compensated 
VAR  filtered  spatial  measurements. 

The  estimate  Qd  of  the  covariance  matrix  Qd  of  the  VAR  filtered  interference  sequence 
can  be  written  as 

Qi  =  ;^  E  .  (2.22) 

Given  3,nd  Qdj  the  estimates  of  {5d/}^^  and  DOA  9d  (if  the  array  manifold  is 

known)  or  ad  (if  the  array  manifold  is  unknown)  can  be  obtained  by  minimizing  the  following 
Weighted  Least  Squares  (WLS)  cost  function, 


1,-1  ^  ^  ^  ^ 

C's  =  E  ~  ^di'Wd/ad)  Qd^  (Scdi  -  bdi'HdiSid)  • 


H  - 


(2.23) 
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Note  that  this  cost  function  is  similar  to  the  one  used  for  the  unstructured  method  in  [1] 
since  the  Fisher  information  matrix  (FIM)  for  dcdi  in  (2.14)  is  also  proportional  to  (see 
[1]  for  more  details). 

Method  1:  To  estimate  6d,  we  must  know  the  array  manifold  a(0d)  as  a  function  of  9d-  The 
method  for  estimating  and  Od  by  using  the  array  manifold  is  referred  to  as  Method 

1.  Without  loss  of  generality,  we  consider  a  uniform  linear  array  (ULA),  where  a(0d)  has  the 
form, 


a(0d) 


1 


J{M-l)^^cosed 


(2.24) 


with  A  being  the  radar  wavelength  and  C  being  the  spacing  between  two  adjacent  sensors. 
Minimizing  C5  in  (2.23)  gives 


bdi  = 


Z  =  0,---,X-1, 


(2.25) 


9d  =  arg  max  ^ 


^<-1  ^  Qidj 


1=0  a^iidiQd^ndiSid 


(2.26) 


Once  9d  is  determined,  is  obtained  with  (2.25)  by  replacing  with  a(0d). 

Method  2:  To  achieve  robustness  against  array  calibration  errors,  we  can  assume  that  Od 
is  completely  unknown.  The  corresponding  method  used  to  estimate  both  and  a^  is 

referred  to  as  Method  2.  Note  that  since  replacing  by  {pbai}fSQ  and  by  a^//?  in 

(2.23),  where  P  is  any  non-zero  complex  scalar,  does  not  change  C5,  and  a^  can  only 

be  determined  up  to  an  unknown  multiplicative  complex  constant.  This  unknown  complex 
scalar,  similar  to  the  unknown  gain  and  initial  phase  of  a  radar  system,  does  not  affect  most 
practical  applications  such  as  ATR.  C5  (see  below). 

Given  minimizing  C5  in  (2.23)  with  respect  to  a^  yields 


.1=0  J  J=0 


(2.27) 
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where  (•)*  denotes  the  complex  conjugate.  Given  a^,  can  be  estimated  by  using 

(2.25).  Hence  we  can  cyclically  iterate  (2.25)  and  (2.27)  to  obtain  the  estimates  of 
and  a^. 

To  start  the  iteration,  we  must  have  an  initial  estimate  of  either  a^  or  {bdi}i'SQ.  Our 
initialization  approach  is  obtained  by  using  the  Singular  Value  Decomposition  (SVD)  [8]. 
Rewriting  C5  in  (2.23),  we  have 


{^di  ^di  —  bdi^  Qdi  (^di  ^di  —  bdi^^  ,  (2.28) 

where 

Qdl  =  /  =  0,  •••,//—  1. 

(2.29) 

To  place  the  most  weight  on  the  term  that  is  associated  with  the  largest  signal  energy, 

we  choose  where  Iq  is  selected  such  that 

is  the  largest  among 

ii::-  - 

a<ii  =  ^d'^di^di, 

(2.30) 

and 

II 

(2.31) 

We  minimize  the  following  approximate  cost  function 

L-1 

0*7  =  53  {^<u  —  bdi^)^  {ocdi  —  bdi^d) , 

(2.32) 

1=0 

which  is  equivalent  to 

Cs  =  1  Ad  -  adbJl^ , 

(2.33) 

where  1|  •  || 

F  denotes  the  FVobenius  norm. 

bd  =  [bdo  •  •  •  bd(L-i)T  > 

(2.34) 

and 

Ad  =  [ttdo  ’  • '  Q!d(L_i)] . 

(2.35) 
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The  C%  in  (2.33)  is  minimized  if 


arfbj  =  u<fi(TdiV^,  (2.36) 

where  u^i  and  v^i  are,  respectively,  the  left  and  right  singular  vectors  associated  with  the 
largest  singular  value  Odi  of  A^.  Then  either  the  initial  estimate  of  b^,  i.e., 

bf  =  (2.37) 

or  the  initial  estimate  of  a^,  i.e., 

4°)  =  W^'udi,  (2.38) 

can  be  used  to  initialize  the  alternating  optimization  approach.  We  use  the  one  in  (2.37)  in 
our  numerical  examples.  The  steps  of  Method  2  are  summarized  as  follows: 

Step  (0):  Obtain  the  initial  estimate  bj°^  of  b^  with  (2.37). 

Step  (1):  Update  with  (2.25)  by  replacing  a^  in  (2.25)  with  the  most  recently 

determined  a^. 

Step  (2):  Update  a^  with  (2.27)  by  replacing  in  (2.27)  with  the  most  recently 

determined 

Step  (3):  Iterate  Steps  (1)  and  (2)  until  practical  convergence,  which  is  determined  by 
checking  the  relative  change  ^  of  the  cost  function  in  (2.23)  between  two  consecutive 
iterations. 

We  remark  that  if  the  range  migration  is  negligible,  i.e.,  r  =  0  in  (2.2),  then 
in  (2.8)  do  not  depend  on  1.  Then  Step  (0)  alone  gives  the  solution  that  minimizes  the  C5 
in  (2.23). 

Summary  of  the  Steps  of  Space-Time  Parameter  Estimation 

The  space-time  parameter  estimates  §d  (or  a^),  b<i}£:i  of  multiple  moving  tar¬ 
gets  can  be  obtained  as  follows: 

Step  I.l:  Assume  D  =  1.  Obtain  {Cbd,  9d  (or  a^),  bd}<i=i  from  yi(n). 

Step  1.2:  Assume  D  =  2.  Compute  y2i(n)  with  (2.13)  using  {d»d,  9d  (or  a^),  b(i}d=i  ob¬ 
tained  in  Step  I.l.  Obtain  {0;^,  9d  (or  a^),  bd}d=2  from  y2i(n).  Next,  compute  yii(n)  with 
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(2.13)  using  {a>d,  Qd  (or  a^),  bd}d:=:2  and  then  re-determine  (or  ^),  bd}d=:i  from 

yij(n).  Iterate  the  previous  two  substeps  until  “practical  convergence”  is  achieved  (to  be 
discussed  later  on). 

Step  1.3;  Assume  D  =  3.  Compute  y3/(n)  with  (2.13)  using  {wd,  6d  (or  a^),  bd}d=:i,2  ob¬ 
tained  in  Step  1.2.  Obtain  {a>d,  Od  (or  a^),  bd}d=3  from  yziin).  Next,  compute  yu(n)  with 
(2.13)  using  {ud,  (or  ad),  bd}d=2,3  and  then  re-determine  {ud,  6d  (or  ad),  bd}d=i  from 
yij(n).  Then,  compute  y2i(^)  with  (2.13)  using  {ud,  Od  (or  ad),  bd}d=i,3  and  re-determine 
{^di  Od  (or  ad),  bd}d=2  from  y2i(^)-  Iterate  the  above  three  substeps  until  “practical  con¬ 
vergence”  is  reached. 

Remain  Steps:  Continue  similarly  until  D  is  equal  to  the  desired  or  estimated  number  of 
targets,  which  is  assumed  to  be  known  or  can  be  determined  by  using  the  Generalized  Akaike 
Information  Criterion  (GAIC)  [5,  9].  d  =  !,••'  ,D. 

The  “practical  convergence”  in  the  iterations  of  the  above  algorithm  can  be  determined 
by  checking  the  relative  change  €  of  the  cost  function  Ci  in  (2.12)  between  two  consecutive 
iterations.  The  steps  leading  to  the  last  step  are  needed  to  provide  good  initial  conditions 
for  the  last  step  of  the  algorithm. 

Estimating  the  target  range  features  {oidk,  fdk}k=i  from  bd,  d  —  1, •••,£>,  is  our  next 
concern. 


2.3.2  Target  Range  Feature  Estimation 

Once  the  sequences  {bdi}f=o,  d  =  1,  •  •  • ,  D,  for  all  targets  are  available,  the  range  feature 
estimates,  {odt,  fdk}k^i,  d,  =  1, •••,!),  can  be  obtained  by  minimizing  the  following  cost 
function, 

C9({aji,  MUi)  =  l|bi  -  (2.39) 

A-.  ^ 

where  bdi  is  the  Ith  element  of  vector  bd, 

«d  =  [Odl  •••  <^dKa\  = 
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and 

F</  =  [f<ii  •  • '  ^dKi]  5  d  =  1,"  •  ,D, 

with 

fdfc  =  [l  ,  d  =  1, •  •  • , ^,  k  =  l,-',Kd. 

The  estimates  {ddit,  fdk}k^i  of  {odfe,  fdk}k^i  can  be  obtained  by  using  the  RELAX  algorithm 
(see  [5]  for  more  details),  which  has  a  similar  structure  as  the  approach  used  for  space-time 
parameter  estimation. 

We  remark  that  our  multiple  moving  target  feature  extraction  algorithm  above  may  have 
used  more  unknowns  than  necessary  at  certain  steps.  We  choose  to  do  so  to  simplify  and 
speed  up  the  algorithm.  For  example,  to  use  the  relaxation-based  optimization  algorithm,  we 
could  estimate  both  the  space-time  and  the  range  parameters  of  the  dth  target  and  subtract 
out  the  dth  target  based  on  the  parameter  estimates  {u>d,  (or  a^),  {ddfc,  fdk}k=i}d=i 
the  iteration  steps.  However,  since  estimating  the  range  parameters  {6;^^,  fdk}k=i  can  be 
computationally  demanding,  we  choose  to  separate  the  range  parameter  estimation  from 
the  space-time  parameter  estimation.  Our  numerical  results  have  shown  little  accuracy 
degradation  but  reduced  computations,  especially  for  large  Kd,  as  a  result  of  the  separate 
space-time  and  range  parameter  estimation. 

2.4  Numerical  Examples 

We  present  several  numerical  examples  to  illustrate  the  performance  of  our  proposed 
algorithm.  In  the  following  examples,  we  assume  that  the  array  is  a  ULA  with  M  =  8;  the 
interelement  distance  between  two  antennas  is  ^  =  A/2;  the  number  of  pulses  in  a  CPI  is 
N  -  16;  the  phase  history  sample  number  is  L  =  16,  i.e.,  an  LRR  range  segment  contains 
16  HRR  range  bins.  Consider  two  targets  {D  =  2)  with  DOAs  6i  =  30°  and  62  =  150°, 
and  Doppler  frequencies  Ui  =  0.2-k  and  0J2  =  0.47r  with  r  =  0.01.  Each  of  the  two  targets 
consists  of  two  closely  spaced  scatterers  {Ki  =  K2  =  2)  with  parameters  an  =  1,  ai2  =  1, 
a2i  =  1,  022  =  1,  fn  =  0.1,  /12  =  0.1  -I-  1/2L,  /21  =  0.3,  and  /22  =  0.3  +  1/2L,  where  the 
subscript  ij  means  the  jth  scatterer  of  the  Rh  target.  The  VAR  filter  order  is  P  =  2.  (No 
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obvious  performance  improvement  is  obtained  by  higher  orders.)  The  number  of  secondary 
range  bins  is  S'  =  4.  We  set  e  =  ^  =  10“^  to  determine  the  “practical  convergence”  in  the 
simulations.  The  mean-squared  errors  (MSEs)  of  the  various  estimates  are  obtained  from 
100  Monte  Carlo  trials. 

We  simulate  the  ground  clutter  as  a  temporally  and  spatially  correlated  Gaussian  random 
process  [10].  The  Clutter-to-Noise  Ratio  (CNR),  defined  as  the  ratio  of  the  clutter  variance 
to  the  noise  variance,  is  set  to  be  CNR=  40  dB.  A  jamming  signal,  which  is  a  zero-mean 
temporally  white  Gaussian  random  process,  also  exists.  The  Jammer- to-Noise  Ratio  (JNR), 
which  is  the  ratio  of  the  jammer’s  temporal  variance  to  the  noise  variance,  is  chosen  as 
JNR=  25  dB  and  the  jamming  impinges  from  Oj  =  120°.  When  array  calibration  errors  exist, 
the  errors  for  different  elements  are  assumed  to  be  independent  and  identically  distributed 
complex  Gaussian  random  variables.  In  our  simulations,  a  complex  Gaussian  random  vector 
with  zero-mean  and  covariance  matrix  0.041  is  added  to  the  array  manifold  to  simulate  array 
calibration  errors,  which  implies  that  the  variance  of  the  calibration  error  for  each  element 
is  0.04. 

We  first  present  an  example  of  no  array  calibration  errors.  Note  that  the  complex  am¬ 
plitude  estimates  {adk}k=i  of  obtained  via  Method  2  are  all  scaled  by  a  common 

unknown  complex  scalar.  To  calculate  their  best  possible  MSEs,  we  scale  them  to  minimize 

C,0  =  lloio  -  Poc£,  (2.40) 


where  ajo  is  the  true  value  of  aa-  Minimizing  (2.40)  with  respect  to  P  yields: 


(8  = 


-  H 

OLd  ttdO 


(2.41) 


Note  that  this  scaling  scheme  is  only  used  to  illustrate  the  complex  amplitude  estimate 
performance;  it  is  not  a  necessary  step  in  a  practical  application  including  ATR  since  only 
the  relative  amplitudes  are  of  interest.  For  comparison  purposes,  the  MSEs  of  the  estimates 
{ddjk}^!  obtained  via  Method  1  are  presented  both  with  and  without  a  scaling  scheme 
similar  to  (2.40). 
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Figures  2.2(a)-2.2(d)  show  the  MSEs  of  the  estimates  of  the  Doppler  frequency,  target 
DOA,  complex  amplitude,  and  range  frequency  of  the  first  scatterer  of  the  first  target  as 
a  function  of  Signal-to-Noise  Ratio  (SNR),  which  is  defined  as  the  ratio  of  |q:uP  to  the 
noise  variance,  and  compare  them  with  the  corresponding  Cramer-Rao  bounds  (CRBs)  (see 
Appendix  B  for  the  CRB  derivation).  (The  results  for  the  other  scatterer  and  the  other 
target  are  similar.)  Due  to  using  the  scaling  scheme  in  (2.40),  the  MSEs  of  {adk}k=i  be 
better  than  the  CRBs  (which  do  not  account  for  such  a  scaling).  We  note  that  as  the  SNR 
increases,  the  MSEs  can  approach  the  corresponding  CRBs,  which  indicates  that  the  clutter 
suppression  scheme  works  well  and  the  parameter  estimation  algorithm  is  highly  accurate. 

Next,  we  present  an  example  when  the  array  calibration  errors  exist.  All  other  parameters 
are  kept  the  same  as  for  Figure  2.2.  Figures  2.3(a)-2.3(d)  show  the  MSEs  of  the  estimates 
of  the  Doppler  frequency,  target  DOA,  complex  amplitude,  and  range  frequency  of  the  first 
scatterer  of  the  first  target,  as  a  function  of  SNR.  Note  that  the  MSEs  of  the  complex 
amplitude  and  range  frequency  estimates  obtained  via  Method  2  are  close  to  the  CRBs  as 
the  SNR  increases,  while  the  MSEs  of  the  complex  amplitude  estimates  for  Method  1  fail  to 
follow  the  CRBs  if  the  scaling  scheme  is  not  used,  though  the  range  frequency  is  still  well 
estimated  via  Method  1.  Hence  from  Figure  2.3,  we  note  that  both  Methods  1  and  2  are 
robust  against  calibration  errors  as  far  as  the  Doppler  frequency,  relative  complex  amplitude, 
and  range  frequency  of  the  scatterer  are  concerned. 

rigid-body  moving  L02  =  0.47r,  DOAs,  and  {bk,  9kt  ^k}k=i,2  primary  data 
We  remark  that  the  above  simulation  results  show  that  Methods  1  and  2  provide  similar 
performances  for  target  relative  complex  amplitude  and  range  frequency  estimation.  Method 
2  avoids  the  1-D  search  over  the  DOA  space  and  usually  requires  only  a  few  (3  6)  iterations. 

To  estimate  the  target  phase  history  sequence  {bu}f~Q  from  the  spatial  signature  estimates 
Method  2  needs  only  about  10%  ~  30%  of  the  amount  of  computations  measured 
in  MATLAB  flops  required  by  Method  1.  (This  difference  is  mainly  due  to  the  fact  that  the 
latter  requires  a  1-D  search  over  the  DOA  space.)  Hence  if  the  array  calibration  errors  are 
known  to  be  significant  enough  to  result  in  a  useless  target  DOA  estimate  or  if  the  target 
DOA  is  not  of  interest.  Method  2  is  preferred  over  Method  1. 
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2.5  Conclusions 


We  have  presented  a  robust  and  accurate  method  for  the  clutter  suppression  and  feature 
extraction  of  multiple  moving  targets  for  airborne  HRR  phased  array  radar.  To  avoid  the 
range  migration  problems  that  occur  in  HRR  radar  data,  we  divided  the  HRR  range  profiles 
into  LRR  segments.  We  have  shown  how  to  use  the  VAR  filtering  technique  to  suppress  the 
ground  clutter  and  use  the  relaxation-based  algorithm  to  extract  the  features  of  multiple 
moving  targets.  The  multiple  moving  target  feature  extraction  problem  is  reduced  to  the 
feature  extractions  of  a  single  moving  target  in  a  relaxation-based  iteration  step.  For  each 
target  and  in  each  iteration,  the  target  phase  history  sequence  and  DOA  (or  the  unknown 
array  manifold)  are  estimated  from  the  spatial  signature  vectors  by  minimizing  a  Weighted 
Least  Squares  (WLS)  cost  function.  Numerical  results  have  demonstrated  that  our  multiple 
moving  target  feature  extraction  algorithm  performs  well  in  the  presence  of  strong  interfer¬ 
ence  including  clutter,  noise,  and  jammer  and  is  robust  against  array  calibration  errors. 

Appendix  A;  Data  Model  Derivation 

We  sketch  below  the  derivation  of  the  data  model  used  in  (2.3).  We  first  establish  the 
data  model  for  the  case  of  a  single  antenna,  then  extend  it  to  the  case  of  the  phased  array 
radar. 

Assume  that  the  radar  transmits  a  group  of  chirp  pulses  with  the  pulse  width  To  and  the 
pulse  repetition  interval  (PRI)  T.  A  normalized  chirp  signal  has  the  form 

s{t)  =  \t\  <  To/2,  (2.42) 

where  fo  denotes  the  carrier  frequency  and  2'y  denotes  the  frequency  modulation  rate.  First, 
we  assume  that  there  are  D  targets  with  different  DOAs,  radial  velocities,  and  RCSs.  At 
time  t,  the  range  of  the  dth  target  is  Rd{t)  =  Rd  +  Vdt,  d  =  1,  •  •  • ,  D,  where  Rd  and  Vd 
denote  the  range  location  and  radial  velocity  of  the  dth  target,  respectively,  when  the  first 
pulse  was  transmitted. 
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Let  t  =  t  —  nT  denote  the  fast  time,  where  n  is  the  pulse  number.  Then  the  received 
signal  from  the  dth  target  is: 

rd{t)  =  a^exp  27r/o  (t-To-  At^  -  +  7  -  tq  -  At^  -  | , 

\i\<  To/2,  n  =  (2.43) 

where  is  determined  by  the  RCS  of  the  dth  target,  tq  =  2Rq/c,  Atj  =  2{Rd  —  Ro)/c  with 
Ro  denoting  a  reference  range  (possibly  corresponding  to  the  center  of  the  target),  c  is  the 
speed  of  light,  and  N  is  the  total  number  of  pulses  transmitted  in  a  CPI.  By  using  s{i  —  tq), 
with  \i  —  tqI  <  To/2,  as  the  reference  signal,  the  dechirped  signal  has  the  form: 


Xd{t,n)  =  rd{t)s*{t- Tq) 

=  Orfexp  -j27r/o  ^ATd  + 

r  /  -  2' 

•  exp  -jj  [  2t  -  2ro  -  Ara - 


Ara  + 


(2.44) 


Since  V^/c  and  Ar^  tend  to  be  very  small  in  practice,  when  the  dwell  time,  {N  —  1)T,  is 
short,  we  can  approximately  express  (2.44)  as: 


Xd{i,n)  fa  a^exp  -j 


AnYd  47roVd 
A  c 


2'yAT^  ij  exp  uTnj 


f  A'iVdT-  ■ 
exp  I  — y - tn 


n  =  0, !,•••, AT— 1, 


(2.45) 


where  A  denotes  the  radar  wavelength,  and  ad  is  did  scaled  by  a  constant  phase.  Note  that, 
in  (2.45)  above,  the  first  exponential  term,  a  linear  function  of  i,  corresponds  to  the  phase 
change  of  the  signal  due  to  the  dth  target  within  a  chirp  pulse,  which  is  caused  by  the  relative 
range  and  target  velocity;  the  second  one  accounts  for  the  phase  shift  (Doppler  shift)  between 
pulses,  which  is  due  to  the  radial  velocity  of  the  dth  target;  finally,  the  last  term  represents 
the  accumulated  phase  shift  from  profile  to  profile,  i.e.,  range  migration. 

We  assume  that  2  (J^max  ~  -Rmin)  /c  ^  Tq,  where  Rmax  and  Rmm  denote  the  maximum 
and  minimum  ranges  between  the  radar  and  target  scatterers,  respectively.  Let  Tg  denote 
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the  sampling  period.  Then  we  can  express  the  sampled  dechirped  signal  as 

x/(n)  =  'Z  n  =  0,  •  •  • ,  iV  -  1,  i  =  0,  •  •  • ,  L  -  1,  (2.46) 

«i=l 

where  L  denotes  the  number  of  data  samples  due  to  each  pulse, 


and 


(2.47) 

\  A  TTC  TT  / 

Ud  =  2irVdT 

(2.48) 

ijVdTTs 

Vd  - 

(2.49) 

Since  both  vj,  and  a>d  depend  on  the  relative  speed  between  the  radar  and  the  target,  if  cOd 
is  known,  so  is  Vd,  and  vice  versa.  Defining  Vd  =  rud,  we  have: 

iTs 


T-  —  = 


Wd  tt/o  -  7To  ’ 


(2.50) 


which  is  independent  of  the  target  velocity. 

For  an  HRR  phased  array  radar  with  the  array  manifold  a(0d),  (2.46)  becomes 

X/  (n)  =  E  a(0d) ,  n  =  0,---,N-l,  Z  =  0,---,L-1.  (2.51) 

d=l 

Here,  we  assume  D  rigid-body  targets  with  Kd  scatterers  for  the  dth  target.  We  also  assume 
that  the  scatterers  of  the  same  target  have  the  same  DOA.  Straightforwardly,  (2.51)  has  the 
form 

Xi(n)  =  E  ( E  e^"‘'"'e^‘^‘'"a(0d),  n  =  0,  •  •  • ,  AT  -  1,  1  =  0,  •  •  • ,  L  -  1,  (2.52) 

d=l  \A:=1  / 

where  fdk,  ^d}  ^-re  the  parameters  of  the  A:th  scatterer  of  the  dth  target.  Finally, 

when  the  interference  is  included,  we  obtain  the  data  model  in  (2.1). 


Appendix  B:  Derivation  of  the  CRBs 

The  CRB  matrix  corresponding  to  the  data  model  in  (2.1)  is  derived  for  the  general  case 
in  which  the  interference  term  ej(n)  includes  clutter,  jamming,  and  noise. 
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Rewrite  (2.1)  as, 


X  =  I]  (lladfcWrffc)  a^  +  e  G  (2.53) 

d=i  \fc=i  J 

where  ®  and  0  denote  the  Kronecker  and  Hadamard  matrix  products,  respectively. 


=  [l  , 

(2.54) 

=  [l  , 

(2.55) 

and 

‘^d=[fJ(0)fJ(l)  •••  fJiL-lf, 

(2.56) 

with 

fd{l)  =  [l  ,  1  =  0,  •  •  • ,  L  -  1. 

(2.57) 

Let  Qjv,  Qc,  and  Qj  be,  respectively,  the  covariance  matrices  of  the  noise, 

and  jamming,  which  are  independent  of  each  other.  Then 

ground  clutter. 

Q  =  £^{ee^}  =  +  Qc  +  Qj, 

(2.58) 

where 

QaT  =  CTn^MNL, 

(2.59) 

Qc  =  ®  Qcj 

(2.60) 

and 

Qj  =  cTj^NL  ®  (a^aj ), 

(2.61) 

with  ajf,  Oq,  and  Oj  being,  respectively,  the  variances  of  the  noise,  clutter,  and  jamming, 
1l  being  the  identity  matrix  of  dimension  L,  Qc  being  as  given  in  [10],  and  aj  denoting  the 
jammer  steering  vector  which  has  the  same  form  as  a(0d)  in  (2.24)  except  that  dd  in  (2.24) 
is  replaced  by  the  jammer  DOA  dj. 


According  to  the  extended  Slepian-Bangs’  formula,  the  ijih  element  of  the  Fisher  infor¬ 
mation  matrix  (FIM)  has  the  form  [11,  12] 


{FIM),  =  Tv(Q-|5Q-‘||)+2Re[g)''Q 


dVi  ^  W)  ’ 


(2.62) 


where 

r  1^ 

Tj  =  Re^(a)  Im^(a)  ct>  0  >  (2.63) 

with  a,  f,  oj,  and  0  being  the  vectors  consisting  of  the  complex  amplitudes,  range  frequencies, 
Doppler  frequencies,  and  DOAs  of  the  targets,  respectively;  denotes  the  derivative  of  x 

with  respect  to  the  ith  parameter  of  ry;  Re(x)  is  the  real  part  of  x  and  Im(x)  is  the  imaginary 
part  of  X.  Note  that  the  FIM  is  block  diagonal  since  the  parameters  in  Q  are  independent 
of  those  in  /x  and  vice  versa.  Hence,  the  CRB  matrix  for  the  target  features  and  the  motion 
parameters  can  be  calculated  from  the  second  term  on  the  right  side  of  (2.62).  Let 


(2.64) 


(2.65) 

(2.66) 

(2.67) 


and 

where  =  [1 
in  degrees.  Let 


gd  =  Cd  [(Wrffc  0  O  Wrf]  0  (dM  O  ad),  (2.68) 

1]^  e  and  Cd  =  -i(7r2/18O)sin(0d7r/18O)  with  9d  being  measured 
G  =  Gf  G(j  Gg  I  )  (2.69) 
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where 


—  Su  ■  ■  ■  Siki  ' ' '  Sdi  ■  ■  ■  z'dKd 


Gf  = 


Ga,  =  jG,,, 

Su  ■  ■  ’  e(ki  ' '  *  Sdi  ' '  ■  Edkd 


and 


Go)  — 


Gb  = 


gr  •••  gs 


gi  •••  Eh 


Then  the  CRB  matrix  for  the  parameter  vector  rj  is  given  by 

CRB(77)  =  [2Re(G^Q-^G)]"^ . 


(2.70) 
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Figure  2.1:  HRR  range  profiles  are  divided  into  LRR  segments  with  each  segment  containing 
L  HRR  range  bins 
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(a) 


(b) 


(c)  (d) 

Figure  2.2;  Comparison  of  MSEs  with  CRBs  as  a  function  of  SNR  for  (a)  Doppler  frequency, 
(b)  DOA,  (c)  amplitude,  and  (d)  range  frequency  of  the  first  scatterer  of  the  first  target  with 
CNR=  40  dB,  JNR=  25  dB,  and  no  array  calibration  error. 
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(a) 


(b) 


(c) 


(d) 


Figure  2.3:  Comparison  of  MSEs  with  CRBs  as  a  function  of  SNR  for  (a)  Doppler  frequency, 
(b)  DOA,  (c)  amplitude,  and  (d)  range  frequency  of  the  first  scatterer  of  the  first  target  with 
CNR=  40  dB,  JNR=  25  dB,  and  array  calibration  error  (covariance  matrix  =  0.041)  . 


27 


3.  Nonparametric  Spectral  Analysis  of  Gapped  Data 
via  an  Adaptive  Filtering  Approach 


3.1  Introduction 

Spectral  analysis  of  incomplete  data  is  needed  in  many  applications  including  astrono¬ 
my,  underwater  acoustics,  radar,  and  communications.  Conventional  fast  Fourier  transform 
(FFT)-based  approaches  have  been  widely  used  for  spectral  analysis  tasks  due  to  their  high 
computational  efficiency.  However,  they  suffer  from  low  resolution  and  poor  accuracy.  Many 
other  spectral  analysis  methods  [1,  2,  3]  with  enhanced  performance  have  been  proposed 
for  spectral  analysis  of  data  sequences.  These  methods  include  parametric  (e.g.  [4]),  non¬ 
parametric  (e.g.  [5,  6,  7]),  and  semi-parametric  (e.g.  [8])  approaches.  In  general,  the 

nonparametric  approaches  are  less  sensitive  to  data  modeling  errors  than  their  parametric 
counterparts.  Most,  if  not  all  of  the  cited  methods,  however,  are  devised  for  sequences  of 
contiguous  measurements  without  any  missing  data  samples  or  gaps. 

The  gapped  (or  missing)  data  problem  usually  arises  when  contiguous  data  measurements 
for  a  long  time  are  hard  to  obtain  or  the  measurements  during  some  intervals  are  not  useful 
due  to  strong  interference  or  jamming  and  must  be  discarded.  For  instance,  in  astronomy 
data  measurements  are  available  only  in  the  form  of  groups  of  samples  separated  by  rather 
long  intervals  for  which  no  reliable  measurements  can  be  taken  (see,  e.g.,  [9, 10, 11, 12]  and  the 
references  therein).  In  astronomical  applications,  special  attention  has  been  paid  to  detecting 
the  presence  of  one  or  more  periodic  signals  from  incomplete  data  measurements  by  using 
techniques  such  as  the  periodogram  [13],  discrete  Fourier  transform  (DFT)  interpolation  [14], 
the  CLEAN  deconvolution  [11],  and  reconstruction  of  evenly  sampled  data  by  making  certain 
assumptions  on  the  data  sequence.  In  radar  signal  processing,  azimuth  angle  diversity  data 
fusion  (see,  e.g.,  [15])  and  data  measurements  in  several  discontinuous  frequency  bands  also 
require  spectral  estimation  for  gapped  or  incomplete  data  sequences. 

In  this  paper,  we  present  an  algorithm  for  nonparametric  complex  spectral  analysis  of 
gapped  data  via  an  adaptive  FIR  filtering  approach,  referred  to  as  the  Gapped-data  Ampli- 
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tilde  and  Phase  Estimation  (GAPES)  algorithm.  The  incomplete  data  sequence  may  contain 
gaps  of  various  sizes.  The  GAPES  algorithm  iterates  the  following  two  steps:  (1)  estimat¬ 
ing  the  adaptive  FIR  filter  and  the  corresponding  complex  spectrum  via  APES  (Amplitude 
and  Phase  Estimation)  [5,  6,  7],  a  nonparametric  adaptive  FIR  filtering  spectral  estimation 
approach,  and  (2)  filling  in  the  gaps  via  an  APES  least  squares  fitting  criterion.  The  initial 
condition  for  the  iteration  is  obtained  from  the  available  data  segments  via  APES. 

The  remainder  of  this  paper  is  organized  as  follows.  In  Section  2,  we  formulate  the 
problem  of  interest.  The  GAPES  algorithm  for  complex  spectral  analysis  of  gapped  data  is 
described  in  detail  in  Section  3.  In  Section  4,  numerical  examples  are  presented  to  illustrate 
the  performance  of  the  proposed  algorithm.  Finally,  Section  5  contains  our  conclusions. 

3.2  Problem  Formulation 

We  consider  the  complex  spectral  analysis  of  one-dimensional  discrete-time  evenly  sam¬ 
pled  data  sequences  with  gaps  of  various  sizes.  Let  {2:(n)}nro^  denote  a  data  sequence 
without  gaps  or  missing  data,  where  N  denotes  the  data  length.  The  spectral  analysis  of 
{a;(n)}  essentially  amounts  to  decomposing  x{n),  at  each  frequency  to  of  interest  as: 

a:(n)  =  a(a;)e^‘^" -f  ea,(n),  n  =  0,  •  •  • ,  AT  —  1,  a;G[0, 27r),  (3.1) 

where  ea,(n)  denotes  the  unmodeled  noise  and  interference  at  frequency  u,  and  a{u})  is  the 
value  of  the  complex  spectrum  of  {a:(n)}  at  ca.  Note  that  in  some  applications  only  the 
amplitude  (or  power)  spectrum  |a(6i;)|  is  of  interest,  whereas  in  others  the  phase  of  a{u)  is 
also  required. 

For  the  complete  data  sequence  {x(n)},  a  simple  method  of  estimating  the  complex 
spectrum  a(a;)  is  to  use  the  computationally  efficient  FFT.  The  FFT  of  {x(n)}  at  frequency 
u  normalized  by  a  factor  N  gives  the  FFT  spectral  estimate  Q:(a;)  of  a(iv)  (note  that  padding 
with  zeros  when  performing  the  FFT  may  be  necessary  to  obtain  a  smooth  complex  spectrum 
[1,  2,  3]).  However,  the  FFT  approach  is  known  to  suffer  from  sidelobe  artifacts  and  poor 
accuracy.  Many  types  of  windows  can  be  applied  to  the  data  sequence  prior  to  performing 
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FFT  to  reduce  the  sidelobe  effects  at  the  cost  of  worsening  the  spectral  resolution. 

When  some  segments  of  the  data  sequence  {x{n)^  are  not  available,  we  are  led  to  the 
gapped  data  problem.  Let 

X  =  [x(0)  •  •  •  x{N  —  1)]^ 

^  [xf  ...  x?f  (3.2) 

be  the  full  data  vector,  where  (.)^  denotes  the  transpose  and  Xi  . . .  xp  are  nonoverlapping 
subvectors  of  x  with  lengths  Ni,  TVp  and  T,p=iNp  =  N.  A  gapped  data  vector  Xg  is 
defined  similarly  to  x  above  except  that  Xp  for  p  =  2,  4,  . .  • ,  P  —  1  (P  is  assumed  to  be  an 
odd  number)  are  missing. 

The  FFT-based  spectral  analysis  of  gapped-data  sequences  by  filling  in  the  gaps  with 
zeros  may  well  yield  an  estimated  spectrum  with  many  artifacts,  such  as  false  peaks  (which 
in  RADAR  applications  would  correspond  to  ghost  targets).  This  phenomenon  is  illustrated 
in  Figure  3.1  in  which  the  simulated  true  spectrum  (Figure  3.1(a))  consists  of  two  spectral 
lines  located  at  fi  =  a;i/27r  =  0.12  Hz  and  /2  =  U2/2'it  =  0.14  Hz  with  amplitudes  ori  =  1 
and  012  =  0.2,  respectively.  The  total  number  of  samples  is  A  =  64  with  the  samples  19 
through  40  missing  and  hence  set  to  zero  in  the  FFT  (in  the  previous  notation  this  scenario 
corresponds  to  Ni  =  18,^2  =  22,  and  A3  =  24).  The  available  data  samples  are  corrupted 
by  an  additive  white  Gaussian  noise  sequence  with  zero-mean  and  variance  0.01.  A  Taylor 
window  with  shape  parameter  5  and  -35  dB  sidelobe  level  will  be  used  throughout  this  paper 
for  the  windowed  FFTs.  Although  usually  we  are  interested  in  the  complex  spectrum,  in  the 
figures  we  only  present  the  moduli  of  the  complex  spectral  estimates  for  convenience.  Figure 
3.1(b)  shows  the  magnitude  of  the  windowed  FFT  spectrum  of  the  gapped  data  sequence 
{xg(n)}  with  a  zero-padding  length  equal  to  4 A.  Owing  to  replacing  the  missing  data  with 
zeros,  the  spectrum  in  Figure  3.1(b)  contains  significant  sidelobe  artifacts,  which  are  due  to 
leakage  from  the  dominant  peak.  The  peak  due  to  the  second  spectral  line  at  /2  =  0.14  Hz 
is  completely  buried  in  the  sidelobe  effects  of  the  dominant  peak.  For  comparison.  Figures 
3.1(c)  and  (d),  respectively,  show  the  moduli  of  the  windowed  FFT  spectra  of  {x(n)}  and 
{x3(n)},  the  longer  non-zero  data  segment  in  {xg{n)}.  Both  sequences  are  zero-padded  to  4 
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times  their  original  lengths.  The  two  spectral  lines  are  not  resolved  in  the  windowed  FFT 
spectrum  in  Figure  3.1(d)  due  to  the  short  length  of  {a:3(n)}.  Furthermore,  even  though 
the  difference  between  the  two  frequencies  fi  and  /2  is  larger  than  the  FFT  resolution  limit 
of  1/N  (i.e.  |/i  -  /2I  =  0.02  >  1/A^  0.016),  they  cannot  be  clearly  resolved  even  when 

the  complete  data  sequence  is  assumed  to  be  available  since  windowing  worsens  the  FFT 
resolution  (see  Figure  3.1(c)). 

Hence,  as  expected,  the  analysis  of  individual  contiguous-data  segments  in  a  gapped 
sequence  leads  to  poor  spectral  resolution.  Also,  data  fusion  via  the  direct  FFT  of  the 
sequence  with  the  gaps  filled  with  zeros  yields  spectra  with  significant  artifacts,  such  as 
ghost  peaks,  which  are  difficult  to  interpret.  In  what  follows  we  show  how  to  estimate 
q:(cc;)  from  {a:p(n)}  via  an  enhanced  nonparametric  adaptive  filtering  approach  referred  to  as 
GAPES  (an  acronym  that  was  explained  in  Section  1). 

3.3  Spectral  Analysis  of  Gapped  Data  via  GAPES 

We  first  present  the  GAPES  algorithm  for  a  data  sequence  with  one  gap.  We  then  extend 
it  in  a  straightforward  manner  to  the  case  of  multiple  gaps  of  various  lengths. 

3.3.1  Data  Sequence  with  One  Gap 

When  only  one  gap  exists  in  the  measured  data  sequence,  we  have 

xj  =  [x[  x^  xf]^,  (3.3) 

where  the  elements  of  X2  are  considered  to  be  unknown  parameters  in  what  follows.  Let 
h(a;)  denote  the  impulse  response  of  an  M-tap  FIR  filter, 

h(a;)  =  [hi(u;)  /i2(a;)  •••  hM{u})]'^,  (3.4) 

and  let 

^gil)  =  [xg{l)  Xg{l  +  1)  •••  Xg{l  +  M  -1)]'^,  i  =  0,---,L-l,  (3.5) 
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denote  the  forward  overlapping  vectors  constructed  from  the  sequence  where  L  = 

N  —  M  +  1  and  Xg{n)  denotes  the  nth  element  of  Xg.  Note  that  Xg{l)  is  a  function  of  X2.  We 
would  like  to  choose  h(a;)  such  that  the  filter  output,  h^Xg{l),  resembles  as  much  as  possible 
the  sinusoidal  component  of  frequency  u  in  the  data  sequence: 

h^(c<;)xp(Z)  =  +  reSij(Z),  Z  =  0,  — 1,  (3.6) 

where  the  residual  term  should  be  as  small  as  possible  (in  a  least-squares  sense,  see  below).  At 
the  same  time,  we  want  to  choose  X2  similarly  so  that  the  {reSu,(/)}  are  kept  small  for  any  cj. 
Finally,  we  determine  the  estimated  complex  spectrum  aiu))  from  the  same  considerations. 
Mathematically,  we  propose  to  obtain  the  estimated  spectrum  as  well  as  an  estimate  of  the 
missing  data  by  solving  the  following  minimization  problem: 

{x2,  h(u;),  0(0;)}  =  min  lli^(‘^)xs(Z)  -  o;(a;)e'‘^'|^  , 

I  J  X2,h(w),a(a.)  ^  ' 

subject  to  h^(a;)a(a;)  =  1,  (3.7) 

where  Y,w  denotes  the  summation  over  a  set  of  frequency  values  in  [0,  2it)  (which  will  be 
specified  later),  and 

a(a;)=[l  •••  .  (3.8) 

In  the  above  equations  we  have  used  the  notation  li(u;)  to  emphasize  the  dependence  of 
the  filter  coefficient  vector  on  the  frequency  oj.  Also,  in  these  equations,  (•)^  denotes  the 
conjugate  transpose,  and  the  constraint  h^(a;)a(w)  =  1  is  imposed  to  ensure  that  the  FIR 
filter  passes  the  frequency  to  undistorted. 

The  minimization  in  (3.7)  is  a  complicated  nonlinear  (fourth-order)  optimization  problem. 
Below  we  present  a  cyclic  approach  that  alternates  between  two  basic  steps  of  spectrum 
estimation  and  gap  filling.  Specifically,  one  step  determines  the  adaptive  filter  and  the 
complex  spectrum  by  using  the  most  recent  estimate  X2  of  X2  in  (3.7).  As  shown  in  the  sequel 
this  step  entails  an  application  of  the  APES  algorithm  to  the  interpolated  data  sequence 
X  =  [xf  x^  The  other  step  fills  in  the  gap  using  the  latest  filter  and  spectrum  estimate 
in  (3.7)  to  determine  X2.  These  two  steps  are  presented  in  detail  below. 
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Spectral  Estimation 


Assume  that  an  estimate  X2  of  X2  is  available.  We  obtain  the  interpolated  data  x  by 
replacing  X2  in  (3.3)  with  X2.  The  estimated  vectors  x(i),  Z  =  0, 1,  •  •  • ,  L  —  1,  are  constructed 
from  X  in  the  same  way  as  Xg{l)  are  formed  from  Xg  (see  (3.5)).  The  criterion  in  (3.7)  for 
filter  and  spectrum  determination  can  then  be  expressed  as 


L-l 


{h(a;),  d(w)}  =  min  |h^(^)x(0  -  Q:(w)e^’ 
subject  to  li'^(a;)a(a;)  =  1. 


ul 


(3.9) 


Let  X(u)  be  the  normalized  Fourier  transform  of  the  vector  sequence  {x(Z)} 


x{w)  =  rES(0e"’"'. 

^  1=0 


and  let  R  denote  the  sample  covariance  matrix  of  {x(i)}, 


(3.10) 


^  1=0 


(3.11) 


Then  it  can  readily  be  verified  [7]  that  the  minimization  of  the  quadratic  criterion  in  (3.9) 
with  respect  to  h(w)  and  a{u))  yields  the  filter 

Q~^(a;)a(cj) 


^  ^  a^(a;)Q“^(t<;)a(w)’ 

and  the  following  estimated  complex  spectrum  at  frequency  w; 


(3.12) 


Q;(a;)  =  h^(a;)X(a;), 


(3.13) 


where 


Q(w)  =  R  -  X(a;)X^(a;).  (3.14) 

The  {q:(w)}  above  is  recognized  as  the  APES  spectrum  estimate  obtained  from  the  inter¬ 
polated  data  sequence  x  [5,  6,  7].  It  might  be  thought  that  the  implementation  of  (3.13) 
requires  the  inversion  of  an  M  x  M  matrix  for  each  u).  However,  this  is  not  necessary.  Indeed, 
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by  the  matrix  inversion  lemma: 


l-X(a;)"R  X(w) 


Hence  we  only  need  to  compute  R  ,  which  is  independent  of  oj.  By  letting 


P{u})  =  a(u;)^R  a(a;) 

7(0;)  =  a(u;)"R“^X(a;) 

p(a;)  =  X(a;)^R"^X(a;) 


(3.15) 


(3.16) 

(3.17) 

(3.18) 


we  can  simply  write: 

and  similarly  for  h(a;).  When  applied  to  complete  data  sequences,  APES  achieves  the  best 
performance  for  a  filter  length  of  M  =  N/2  [5,  6].  In  what  follows,  we  will  choose  M  = 
N/2  (unless  specified  otherwise)  to  calculate  the  APES-like  spectrum  estimate  in  (3.13) 
for  the  gapped-data  case.  We  will  evaluate  (3.13)  at  the  discrete  Fourier  transform  (DFT) 
frequencies  cok  =  2Tik/K,  k  =  0,1,  -  •  •  K  —  1  (with  K  >  N). 


Gap  Filling 

Once  {h(wjfc),  6:(wfc)}fro^  were  made  available,  we  can  determine  X2  by  solving  the  fol¬ 
lowing  quadratic  minimization  problem  (cf.  (3.7)): 

(3.20) 

h^(wjk) 
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and 


Zfc  =  a{u)k) 


p3^k 


Using  this  notation  we  can  write  the  criterion  in  (3.20)  as: 

Il2 


K-\ 

E 

fc=0 


Xi 

Hfc  X2  -  Zfc 

X3 


Partition  Hfc  as: 


and  define 


Hfc  =  ^  ^ 

Ni  N2  Nz 


Yk  =  -  AfcXi  -  CfcXs  +  Zfc. 


With  this  notation  the  criterion  becomes, 


K-l 


-  YfclP 


fc=0 


whose  minimizer  with  respect  to  X2  is  well-known  to  be 

■K-l  1-1  K-l 


X2  = 


K-l  1  K-l 

s  ■£  Bfy». 

fc=0  J  fc=0 


(3.22) 


(3.23) 


(3.24) 


(3.25) 


(3.26) 


(3.27) 


The  interpolated  data  sequence  x  is  obtained  by  replacing  X2  in  (3.3)  with  X2. 


Initialization 

To  start  the  algorithm  we  need  an  initial  estimate  of  either  X2  or  the  adaptive  filter  along 
with  the  complex  spectrum.  In  our  initialization  approach  we  first  obtain  an  initial  adaptive 
FIR  filter  h(°)(a;)  and  an  initial  spectral  estimate  a(‘’)(a;)  by  using  the  available  data  Xi  and 
X3. 
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Note  that  both  the  forward  and  backward  data  vectors  are  often  used  to  achieve  better 
spectral  estimation  performance  [5,  6].  When  both  the  forward  and  backward  vectors  are 
used  for  spectral  estimation,  h(a;)  in  (3.12)  and  6;(tu)  in  (3.13)  remain  the  same  except  that 
Q(t(;)  in  (3.14)  is  replaced  by  [7] 

Q(u;)  =  R-X(cj)X^(a;),  (3.28) 

where 


R  = 

(3.29) 

X(a;)  = 

i[x(a,)  X(u.)], 

(3.30) 

X(u,)  = 

^  1=0 

(3.31) 

R  = 

(0, 

^  1=0 

(3.32) 

and  where  x(/),  /  =  0, 1,  •  •  • ,  L  —  1  are  the  backward  overlapping  vectors 

k{l)  =  [x*{N  -l-l)  x*{N-l-2)  •••  x*iN-l-M)f.  (3.33) 


The  filter  length  we  use  in  the  initialization  step  is  the  smallest  integer  larger  than  or 
equal  to  0.75min(7Vi,  N3).  The  h(®^(a;)  is  calculated  as  in  (3.12)  in  which  Q(®)(ti;)  is  obtained 
by  using  (3.28)  with  the  following  modifications: 


1  /^(O)  =(0)  -(0)  ::(0) 

R'”>  =  i(R;  +R,  +R,  +R^ 


(3.34) 


and 

X(")(w)  =  1  [xf>(a.)  X^">(«)  xi^Hu)],  (3.35) 

i;  (0)  £  (0) 

where  and  R^  ,  z  =  1,3,  are  the  estimated  forward  and  backward  sample  covariance 
matrices  obtained  from  Xj,  i  =  1,3;  and  similarly  xj°^(6i;)  and  X^^^o;),  i  =  1,3,  are  the 
normalized  Fourier  transforms  of  the  forward  and  backward  overlapping  vectors.  To  compute 
those  quantities  we  use  their  previously  given  definitions  with  L  replaced  by  Lj  =  ATj  — M  +  1. 
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The  initial  APES  spectrum  is  given  by 

dW(wfc)  =  [^iM  +  ,  (3-36) 

for  k  =  0,1,- •  •  ,K  —  1. 

Summary  of  the  Steps  for  One  Gap 

Step  0.  Initialization:  start  the  algorithm  using  the  initialization  approach  presented  above. 
Step  1.  Gap  filling:  update  the  gap  estimate  by  using  (3.27)  based  on  the  most  recently 
determined  adaptive  filter  and  complex  spectrum  {h(a;fe),  Q:(wfc)}. 

Step  2.  Spectral  estimation:  determine  the  adaptive  filter  and  the  complex  spectrum  with 
(3.12)  and  (3.13),  respectively,  from  the  most  recently  interpolated  data  sequence  {x(n)}. 
Step  3.  Iteration:  repeat  Steps  1  and  2  until  the  relative  change  of  the  cost  function 

2  |h»K)S(i)  -  a(w»)e'“‘'f  (3-37) 

fc=0  l-O 

between  two  consecutive  iterations  is  equal  to  or  less  than  a  preset  threshold  value  c  or  the 
number  of  iterations  reaches  a  preset  number  /max- 

3.3.2  Data  Sequence  with  Multiple  Gaps 

The  extension  of  the  algorithm  presented  in  the  previous  section  to  the  multi-gap  case 
is  straightforward.  The  criterion  (3.7)  that  lies  at  the  basis  of  GAPES  remains  unchanged 
and  so  does  the  spectrum  estimation  step  of  GAPES.  The  gap  filling  step  requires  a  slight 
modification  of  the  partition  in  (3.24)  to  take  into  account  the  multi-gap  structure  of  the 
data,  and  so  does  the  initialization  step. 

In  fact  it  should  be  clear  from  the  previous  discussion  that  GAPES  is  applicable  to  cases 
in  which  the  missing  data  has  any  pattern.  The  only  step  that  may  need  a  more  significant 
modification  in  such  cases  is  the  initialization  step.  In  the  case  of  arbitrary  missing-data 
patterns  we  can  hardly  use  any  special  initialization.  Hence  we  may  have  to  use  the  simplest 
initialization  of  all  possible  initializations:  set  missing  data  to  zero.  If  the  ratio  of  the  number 
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of  available  data  samples  to  that  of  missing  data  is  reasonably  high  the  aforementioned  simple 
initialization  may  be  quite  satisfactory. 

Finally,  we  remark  on  the  fact  that  the  GAPES  algorithm  can  be  readily  extended  to 
criteria  of  the  more  general  form 

X!  Z)  |h^(a;it)xp(Z)  -  Q;(wfc)e'^'f  (3.38) 

k=0  1=0 

where  the  weights  {wjc  >0;  >  0}  can  be  used  to  emphasize  a  frequency  band  in  which 

we  have  a  particular  interest  or  a  subsample  of  the  data  sequence  that  is  thought  of  being 
more  reliable. 

3.4  Numerical  Results 

We  present  several  numerical  examples  to  illustrate  the  performance  of  GAPES  for  the 
complex  spectral  analysis  of  gapped  data.  We  choose  K  =  N  and  M  =  N/2  in  the  iteration 
steps  of  GAPES.  Although  the  frequency  grid  for  &{(v)  can  be  chosen  in  many  ways,  we  use 
the  DFT  grid  cjk  =  2'KklK,  k  =  0,1,  -  •  •  ,K  —  1.  The  noise  added  to  the  available  data 
samples  is  white  Gaussian  with  zero-mean  and  variance  0.01.  In  the  initialization  step  of 
GAPES,  the  filter  length  is  chosen  to  be  the  smallest  integer  larger  than  or  equal  to  3/4  of 
the  length  of  the  shortest  available  data  segment.  We  use  c  =  10“^  to  test  the  convergence 
of  the  GAPES  algorithm  and  set  the  maximum  number  of  iterations  to  Imax  =  10.  Once 
the  iterative  process  is  terminated,  we  obtain  our  final  spectral  estimate  by  evaluating  the 
estimated  spectrum  of  the  last  iteration  on  a  grid  that  is  four  times  finer  than  the  Jf-point 
DFT  grid. 

We  compare  our  GAPES  algorithm  with  the  windowed  FFT  approach  as  well  as  the 
parametric  CLEAN  algorithm  (see,  e.g.,  [10,  11])  that  has  been  proved  to  be  quite  useful 
in  many  applications  including  astronomical  data  analysis,  microwave  imaging,  and  target 
feature  extraction.  When  using  CLEAN,  we  model  the  data  as  a  sum  of  complex  sinusoids. 
CLEAN  first  estimates  the  parameters  of  the  strongest  sinusoid  via  FFT  with  zero-padding, 
subtracts  it  from  the  original  data,  and  then  repeats  this  process  for  the  next  strongest 
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sinusoid  until  the  predetermined  number  of  sinusoids  is  reached  or  the  estimated  amplitude 
of  the  latest  sinusoid  is  small  enough  compared  to  the  strongest  one.  Once  the  sinusoidal 
parameters  are  estimated,  they  are  used  to  simulate  the  data  in  the  gaps.  The  windowed 
FFT  with  zero-padding  is  then  applied  to  the  so-obtained  interpolated  data  sequence  to  form 
the  CLEAN  spectrum. 

For  the  same  example  as  the  one  in  Figure  3.1,  Figure  3.2(a)  shows  the  modulus  of  the 
APES  spectrum  obtained  from  the  complete  data  sequence  {x(n)}.  A  comparison  between 
Figures  3.2(a)  and  3.1(c)  shows  that  APES  outperforms  the  windowed  FFT  approach.  Fig¬ 
ures  3.2(b)  and  (c)  show  the  moduli  of  the  windowed  FFT  spectra  obtained  from  the  data 
interpolated  by  CLEAN  and  GAPES,  respectively.  Note  that  in  Figure  3.2(b)  the  weaker 
spectral  line  at  =  0.14  Hz  is  not  visible.  Figure  3.2(c)  is  comparable  to  Figure  3.1(c) 
which  indicates  that  the  interpolated  data  obtained  by  GAPES  are  close  to  the  true  values 
of  the  unavailable  data.  Figures  3.2(d)  shows  the  modulus  of  the  APES  spectrum  of  the  data 
interpolated  by  GAPES,  which  is  almost  as  good  as  that  in  Figure  3.2(a)  obtained  from  the 
complete  data  sequence. 

Next  consider  a  gapped  data  sequence  with  the  same  length  and  the  same  gap  size  and 
location  as  in  the  previous  example  but  with  a  more  complicated  spectrum.  Specifically,  a 
mbced  spectrum  (shown  in  Figure  3.3(a))  is  considered  in  this  example,  which  comprises  two 
discrete  spectral  lines  located  at  /i  =  0.05  Hz  and  /a  =  0.07  Hz  with  amplitudes  ai  =  l  and 
CKa  =  0.5,  respectively,  and  a  continuous  spectral  component  centered  at  fa  =  0.18  Hz  with 
a  width  of  6  =  0.04  Hz  and  a  constant  magnitude  of  a:  =  0.25.  Figures  3.3(b)  and  (c)  show 
the  moduli  of  the  windowed  FFT  spectra  of  the  gapped  and  complete  data  sequences.  The 
moduli  of  the  windowed  FFT  spectra  obtained  from  the  CLEAN  and  GAPES  interpolated 
data  sequences  are  presented  in  Figures  3.3(d)  and  (e).  Figure  3.3(f)  shows  the  magnitude 
of  the  APES  spectrum  of  the  GAPES  interpolated  data  sequence.  Comparisons  between  the 
various  parts  of  Figure  3.3  lead  to  similar  conclusions  to  those  drawn  from  the  results  shown 
in  Figures  3.1  and  3.2. 

Finally,  we  present  an  example  that  concerns  the  spectral  analysis  of  an  incomplete  data 
sequence  with  two  gaps.  The  total  length  of  the  simulated  sequence  is  N  =  128,  and  the 
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samples  23  through  46  and  76  through  100  are  missing.  Hence  nearly  40%  of  the  data  samples 
are  missing.  Figure  3.4(a)  shows  the  true  spectrum  which  consists  of  four  discrete  spectral 
lines  at  /i  =  0.05  Hz,  ^  =  0.065  Hz,  /a  =  0.26  Hz,  and  =  0.28  Hz  with  amplitudes 
Q-j  =  q;2  =  <^3  =  1  a^nd  0:4  =  0.5,  and  a  continuous  spectral  component  at  /*  =  0.18  Hz  with 
a  width  b  =  0.015  Hz  and  a  constant  magnitude  of  a  =  0.25.  Figures  3.4(b)  and  (c)  show  the 
moduli  of  the  windowed  FFT  spectra  of  the  gapped  and  complete  data  sequences.  Figures 
3.4(d)  and  (e)  display  the  moduli  of  the  windowed  FFT  spectra  of  the  interpolated  data 
sequences  obtained  via  CLEAN  and  GAPES,  respectively.  Figure  3.4(f)  shows  the  modulus 
of  the  APES  spectrum  of  the  interpolated  data  sequence  obtained  via  GAPES.  Note  that 
in  this  case  CLEAN  cannot  remove  all  false  peaks  owing  to  the  serious  leakage  problem  of 
FFT,  which  makes  CLEAN  pick  peaks  at  wrong  locations.  Once  again,  GAPES  appears  to 
provide  the  most  accurate  interpolated  data  sequence  and  spectral  estimate. 

3.5  Conclusions 

We  have  presented  an  algorithm  for  nonparametric  complex  spectral  analysis  of  incom¬ 
plete  data  with  gaps  of  various  lengths  via  an  adaptive  FIR  filtering  approach,  referred  to 
as  the  Gapped-data  Amplitude  and  Phase  Estimation  (GAPES)  algorithm.  The  GAPES 
algorithm  iterates  the  following  two  steps:  (1)  estimating  the  adaptive  filter  and  the  cor¬ 
responding  complex  spectrum  via  APES,  and  (2)  filling  in  the  gaps  by  using  the  APES 
least-squares  fitting  criterion.  We  initialized  the  GAPES  algorithm  by  applying  APES  to 
the  available  data  segments.  Numerical  examples  have  been  used  to  demonstrate  the  effec¬ 
tiveness  of  the  proposed  GAPES  algorithm  and  compare  it  with  the  windowed  FFT  approach 
and  the  parametric  CLEAN  algorithm.  The  results  have  shown  that  GAPES  yields  more 
accurate  interpolated  sequences  and  spectral  estimates  than  the  other  approaches  considered 
in  this  chapter. 
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Figure  3.1:  Moduli  of  the  complex  spectral  estimates  (solid  lines  in  (b)-(d))  compared  with 
the  true  spectrum  (dotted  lines  in  (b)-(d));  (a)  true  spectrum,  (b)  windowed  FFT  spectrum 
of  the  gapped  data  sequence,  (c)  windowed  FFT  spectrum  of  the  complete  data  sequence, 
(d)  windowed  FFT  spectrum  of  the  longer  contiguous-data  segment. 
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(c)  (d) 

Figure  3.2;  Moduli  of  the  complex  spectral  estimates  for  the  same  example  as  in  Figure 
3.1.  (a)  APES  spectrum  of  the  complete  data  sequence,  (b)  windowed  FFT  spectrum  of 
the  interpolated  data  sequence  obtained  via  CLEAN,  (c)  windowed  FFT  spectrum  of  the 
interpolated  data  sequence  obtained  via  GAPES,  (d)  APES  spectrum  of  the  interpolated 
data  sequence  obtained  by  GAPES. 
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(e)  (f) 

Figure  3.3;  Moduli  of  the  complex  spectral  estimates  (solid  lines  in  (b)-(f))  compared  with  the 
true  spectrum  (dotted  lines  in  (b)-(f));  (a)  true  spectrum,  (b)  windowed  FFT  spectrum  of  the 
gapped  data,  (c)  windowed  FFT  spectrum  of  the  complete  data,  (d)  windowed  FFT  spectrum 
of  the  interpolated  data  sequence  obtained  via  CLEAN,  (e)  windowed  FFT  spectrum  of  the 
interpolated  data  sequence  obtained  via  GAPES,  and  (f)  APES  spectrum  of  the  interpolated 
data  sequence  obtained  via  GAPES.  45 
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Figure  3.4:  Moduli  of  the  complex  spectral  estimates  (solid  lines  in  (b)-(f))  compared  with  the 
true  spectrum  (dotted  lines  in  (b)-(f));  (a)  true  spectrum,  (b)  windowed  FFT  spectrum  of  the 
gapped  data,  (c)  windowed  FFT  spectrum  of  the  complete  data,  (d)  windowed  FFT  spectrum 
of  the  interpolated  data  sequence  obtained  via  CLEAN,  (e)  windowed  FFT  spectrum  of  the 
interpolated  data  sequence  obtained  via  GAPES,  and  (f)  APES  spectrum  of  the  interpolated 
data  sequence  obtained  via  GAPES.  46 


4.  Spectral  Estimation  of  Gapped  Data  and  SAR 

Imaging 

with  Angle  Diversity 


4.1  Introduction 

Estimating  the  spectrum  of  a  time  series  is  crucial  in  applications  such  as  radar,  commu¬ 
nications,  underwater  acoustics,  and  astronomical  data  analysis.  Conventional  fast  Fourier 
transform  (FFT)  based  approaches  have  been  widely  used  due  to  their  high  computational 
efficiency.  However,  the  major  drawback  of  the  FFT-based  approaches  is  the  inherent  low 
resolution.  A  large  number  of  methods  have  been  proposed  to  obtain  enhanced  spectral 
estimates  with  higher  resolution.  Of  these,  the  parametric  (e.g.  [1,  2,  3,  4]),  non-parametric 
(e.g.  [5,  6,  7,  8]),  and  semi- parametric  [9]  approaches  have  been  considered  for  target  feature 
extraction  and  radar  imaging.  In  general,  nonparametric  approaches  are  less  sensitive  to 
data  modeling  errors  and  more  robust  than  their  parametric  counterparts.  Many  nonpara¬ 
metric  spectral  estimators  make  use  of  adaptive  (data-dependent)  finite  impulse  response 
(FIR)  filterbanks.  The  APES  (Amplitude  and  Phase  Estimation)  spectral  estimator  [7,  10] 
belongs  to  this  class  of  approaches. 

The  gapped-data  (or  missing-data)  problem  arises  when  continuous  measurements  are 
not  possible,  or  the  measurements  during  certain  periods  are  not  valid  due  to  for  instance 
interference  or  jamming  impinging  on  the  receiver  and  hence  they  must  be  discarded.  In 
radar  signal  processing  (see,  e.g.,  [11,  12]),  the  problem  of  combining  several  sets  of  mea¬ 
surements  made  at  different  azimuth  angle  locations  can  be  posed  as  a  problem  of  spectral 
estimation  from  gapped  data.  Similar  problems  arise  in  data  fusion  via  ultrawide-band  co¬ 
herent  processing  [13].  Furthermore,  in  astronomy,  data  are  often  available  as  groups  of 
samples  with  rather  long  intervals  during  which  no  measurements  can  be  taken  (see,  e.g., 
[14,  15,  16,  17,  18,  19]  and  the  references  therein).  In  astronomy,  special  attention  has  been 
paid  to  detecting  the  presence  of  one  or  more  periodic  signals  from  incomplete  measurements 
using  techniques  such  as  the  periodogram  approach  [15]  and  the  CLEAN  deconvolution  [16]. 
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In  this  paper,  we  consider  using  the  APES  filter  for  the  spectral  estimation  of  gapped  data 
and  synthetic  aperture  radar  (SAR)  imaging  with  angle  diversity.  Unlike  the  FFT  and  the 
windowed  and/or  averaged  FFT  spectra,  the  APES  spectrum  has  good  resolution  properties, 
suffers  from  little  or  no  leakage  effects,  and  has  good  statistical  stability  [7,  20].  The  excellent 
performance  of  APES  in  the  class  of  nonparametric  spectral  analysis  methods  was  one  of 
the  reasons  why  we  chose  to  extend  this  particular  approach  to  the  gapped-data  case.  We 
present  in  the  paper  a  relaxation  based  algorithm,  referred  to  as  the  GAPES  (Gapped-data 
APES)  algorithm,  which  consists  of  two  steps:  (1)  estimating  the  adaptive  filter  and  the 
corresponding  spectrum  via  APES,  and  (2)  filling  in  the  gaps  via  a  least  squares  (LS)  fit. 
For  the  SAR  imaging  with  angle  diversity  data  fusion,  we  perform  the  one-dimensional  (1-D) 
windowed  FFTs  (WFFTs)  in  range,  use  the  GAPES  algorithm  to  interpolate  the  gaps  in 
the  aperture  for  each  range,  apply  the  1-D  inverse  FFTs  (IFFTs)  and  de-window  in  range, 
and  apply  the  two-dimensional  (2-D)  APES  algorithm  to  the  interpolated  matrix  to  obtain 
the  final  2-D  SAR  image. 

The  remainder  of  this  paper  is  organized  as  follows.  In  Section  2,  we  formulate  the 
problem  of  interest.  The  APES  algorithm  for  the  spectral  analysis  of  the  complete  data  is 
introduced  in  Section  3.  The  GAPES  algorithm  for  the  spectral  analysis  of  the  1-D  gapped 
data  is  presented  in  Section  4.  Section  5  presents  a  method  to  use  GAPES  for  2-D  SAR 
imaging  with  angle  diversity  data  fusion.  In  Section  6,  numerical  examples  are  presented 
to  illustrate  the  performance  of  the  proposed  algorithm.  Finally,  Section  7  contains  our 
conclusions. 

4.2  Problem  Formulation 

We  consider  herein  the  spectral  estimation  algorithm  for  1-D  data  sequences  with  multiple 
gaps  of  various  sizes,  and  for  2-D  data  matrices  with  gaps  consisting  of  missing  columns. 

Let  {x(n)}^Jo  denote  a  complete  1-D  discrete-time  sequence  of  length  N.  The  spectral 
analysis  of  {a:(n)}  essentially  amounts  to  decomposing  a:(n),  at  each  frequency  iv  €  [0,  27r) 


48 


of  interest  as 


x{n)  =  Q:(a;)e^‘^"  +  e^in),  n  =  0,  •  •  • ,  N  —  1,  (4.1) 

where  a!(a;)  denotes  the  complex  amplitude  of  a  sinusoid  with  frequency  w,  and  ea,(n)  denotes 
unmodeled  noise  and  interference  at  u)  (in  other  words,  e,j{n)  is  everything  left  in  x{n)  after 
the  sinusoidal  component  Q;(c<;)e^‘^"  has  been  subtracted).  Note  that  the  decomposition  in 
(4.1)  is  frequency  dependent  and  it  should  not  be  confused  with  assuming  that  x{n)  consists 
of  one  sinusoid  in  noise.  In  fact  Equation  (4.1)  does  not  imply  any  modeling  assumptions  at 
all,  which  means  that  the  methods  based  on  (4.1)  are  nonparametric. 

Let  {x(n,  n)},  n  =  0, 1,  •  ■  • ,  iV  -  1,  n  =  0, 1,  •  •  • ,  iV  -  1,  denote  a  2-D  discrete-time  data 
matrix.  For  a  frequency  pair  of  interest,  {a:(n,n)}  is  modeled  as 

x{n,  n)  =  aipj,  fi), 

n  =  0,  •  •  •,iV  -  1,  n  =  0,  •  •  • ,  JV  -  1;  a;,a;  €  [0,27r),  (4.2) 

where  Q;(a»,  u)  denotes  the  complex  amplitude  of  a  sinusoid,  and  n)  denotes  unmodeled 

noise  and  interference  at  frequency  (a;,cD). 

Assume  that  some  segments  of  the  data  sequence  {x{n)}nZQ  are  unavailable  for  reasons 
explained  above.  Let 

X  =  [a;(0)  •  •  •  x{N  —  1)]^ 

^  [xT  xl’  x?f  (4.3) 

be  the  complete  data  vector,  where  xi,  •••xp  are  subvectors  of  x,  whose  lengths  are 
Ni,“ '  ,Np,  respectively,  with  Y^p=iNp  =  N.  Here  (•)^  denotes  the  transpose.  A  gapped- 
data  vector  Xg  is  formed  when  Xp  for  p  =  2,  4,  •  •  • ,  F  —  1  (F  is  always  an  odd  number) 
are  unavailable  and  considered  as  being  unknown.  Let  Xa  and  Xu  denote  the  two  vectors  of 

dimension  Ni  -f  A3  H - h  Np  and  Az  H - f-  Ap_i,  whose  elements  are  the  available  and 

unavailable  data  samples  in  x,  respectively.  Hence  Xa  is  given  whereas  Xu  is  unknown. 

For  the  angle  diversity  data  fusion  problem  encountered  in  SAR  imaging,  the  data  mea¬ 
sured  by  the  radar  are  recorded  intermittently.  Assuming  that  the  radar  measurements  in 
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range  are  complete,  the  gapped  data  due  to  the  angle  diversity  can  be  expressed  as; 


x^{n,n)  =  < 


unknown, 
x{n,  n), 


n  =  0, 1,  •  •  • ,  iV  —  1;  tieQ, 
otherwise. 


(4.4) 


where  Q  denotes  the  set  of  azimuth  look  indices  n  at  which  no  measurements  are  taken. 

Figure  4.1  shows  a  schematic  map  of  a  spotlight  SAR  with  angle  diversity.  For  two 
adjacent  continuous  segments  and  the  gap  between  them,  we  will  use  the  two  continuous 
segments  to  interpolate  the  gap  between  them.  The  motivation  for  this  is  that  in  many  cases 
the  measurements  can  be  non-stationary. 

The  FFT  processing  for  the  gapped  data  can  yield  a  spectrum  with  significant  artifacts, 
as  demonstrated  and  explained  in,  e.g.,  [21].  The  approach  taken  in  this  paper  is  to  fill  in 
the  gaps  corresponding  to  the  missing  data  and  use  the  so-obtained  interpolated  data  for 
spectral  estimation  via  APES. 


4.3  APES  for  Complete  Data 

We  review  the  basic  facts  about  APES  for  complete  data  [7,  10]  that  are  essential  to 
understanding  the  GAPES  algorithm  for  gapped  data  presented  in  the  next  section.  For 
clarity  of  presentation,  we  discuss  only  the  so-called  forward-only  version  of  APES  (for  both 
1-D  and  2-D  cases).  It  can  be  shown  [7]  that  a  slightly  modified  approach,  leading  to  the 
forward-backward  APES,  gives  significantly  better  estimation  accuracy.  For  this  reason,  the 
forward-backward  APES  will  be  used  in  the  numerical  examples  in  Section  6. 

Estimation  of  0(0;)  in  (4.1)  by  the  method  of  least  squares  (LS)  leads  to  the  Fourier 
transform  (FT)  method: 

n=0 

or  its  counterpart  for  power  spectra,  the  periodogram.  As  is  well-known,  the  FT  approach 
suffers  from  leakage,  poor  resolution  and  erratic  statistical  behavior  and  it  is  rarely  used 
without  some  corrections  [22,  23,  24].  Smoothing  is  one  of  the  most  common  corrections 
applied  to  the  FT.  One  form  of  smoothing  (or  averaging)  of  the  FT  that  will  help  us  make 


50 


the  connection  with  APES  is  as  follows.  Partition  the  data  sequence  in  (4.1)  into  M  x  1 
overlapping  vectors  with  the  following  shifted  structure: 

^{l)  =  [x{l)  ...  x{l  +  M-l)],  /  =  (4.6) 


where  L  =  N  —  M  +  1  and  M  is  a  user  parameter  referred  to  as  the  filter  length.  Then 
calculate  the  FT  of  each  of  the  data  snapshots  in  (4.6)  and  average  the  so-obtained  local 
FTs: 


■t  L-l  1  M-1 
^  Z=0  m=0 


L-l 


^  /=0 


1  Af— 1 

-  a:(Z  4- 

m-0 


(4.7) 


To  interpret  (4.7)  in  a  filterbank  framework,  define 


a(w)  =  [l  e>"  ...  . 


(4.8) 


Observe  from  (4.1)  and  (4.6)  that  we  have 

x(/)  =  -t-  §0,(0)  Z  =  0, . . . ,  L  —  1, 


(4.9) 


where  e^{l)  is  formed  from  {eu{n)}  in  the  same  way  as  x(Z)  is  formed  from  {a;(n)}  in  (4.6). 
Let 

h(a;)  =  ^a(a;).  (4.10) 

The  inner  sum  in  (4.7)  is  equal  to  h^(a;)x(Z),  where  (•)^  denotes  the  conjugate  transpose, 
and  hence: 

q:aft(‘^)  =  Y  [h"(a,)5i(I)]  (4.11) 

/=0 

where  (cf.  (4.9)) 

h"(a;)x(Z)  =  a:(u;)eJ‘"'4-res.,  Z  =  0,...,L-1.  (4.12) 

The  averaged  FT  (AFT)  spectral  estimator  in  (4.7)  can  thus  be  interpreted  as  consisting 
of  two  steps:  (1)  use  the  (frequency-dependent)  finite-impulse-response  (FIR)  filter  h((u) 
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to  prefilter  the  data  sequence  (cf.  (4.12)),  and  (2)  obtain  an  estimated  spectrum  from  the 
so-obtained  filtered  sequence  by  the  method  of  LS,  which  obviously  leads  to  (4.11). 

The  filterbank  interpretation  of  AFT  is  useful  since  it  allows  us  to  understand  the  de¬ 
ficiencies  of  this  approach  by  simply  studying  its  associated  filterbank  {h(a;)}.  As  h(a;)  is 
data-independent  (or  non-adaptive)  it  is  no  surprise  that  d;AFT(^^)  suffers  from  leakage  and 
poor  resolution  problems.  A  well-designed  data- dependent  filterbank  should  be  able  to  re¬ 
duce  such  problems  significantly.  Indeed,  such  a  filter  should  be  able  to  achieve  a  much  larger 
signal-to-interference-and-noise  ratio  (SINK)  in  the  filtered  data  (4.12)  than  the  SINK  in  the 
original  sequence  in  (4.1).  If  this  increase  in  SINR  is  significant  enough  to  counterbalance 
the  reduction  in  the  number  of  data  samples  from  N  for  (4.1)  to  L  for  (4.12),  then  the  result 
will  be  a  spectral  estimate  Q:(a;)  with  enhanced  accuracy.  This  is  precisely  the  idea  behind 
APES,  as  explained  next. 

In  the  APES  approach  the  filterbank  {h(a;)}  is  designed  such  that:  (a)  the  filtered 
sequence  is  as  close  to  a  sinusoidal  signal  as  possible  (in  a  LS  sense),  and  (b)  the  complex 
spectrum  a;(a;)  of  x(n)  is  not  distorted  by  filtering.  Mathematically,  we  obtain  h(a;)  along 
with  the  estimate  of  a(a;)  by  minimizing  the  following  LS  criterion: 


min  Y)  |h^(a))x(Z)  -  a{(ja)e'‘^^f ,  subject  to  h^(a;)a(a;)  =  1,  (4.13) 

a(u)),h(w)  I 


where  the  constraint  h^(a;)a(a;)  =  1  is  imposed  to  ensure  that  h(a;)  passes  a  sinusoid  with 
frequency  u  undistorted  (see  requirement  (b)  above).  The  solution  to  the  above  quadratic 
minimization  problem  in  (4.13)  is  readily  obtained  as  [10] 


h(u;) 


Q  ^(a;)a(a;) 
a»(a;)Q-Hw)a(a;)’ 


(4.14) 


and 

Q;(a;)  =  h^(a;)X(a;), 


(4.15) 


where  X(6(;)  is  the  following  normalized  FT: 

^  1=0 


(4,16) 
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and 


Q(w)  =  R-X(a;)X"(a)) 


(4.17) 


with  R  being  the  following  sample  covariance  matrix: 

(4.18) 

^  /=0 

To  avoid  the  inversion  of  an  M  x  M  matrix  for  each  lo,  we  use  the  matrix  inversion  lemma 
(see,  e.g.,  [22])  to  obtain: 


l-X^^(a;)R  X(a;) 

(4.19) 

-  -1/2  £  -1 

Let  R  denote  the  Cholesky  factor  of  R  ,  and  let 

-1/2 

a(a;)  =  R  a(a;) 

(4.20) 

V.  -1/2  _ 

X(a;)  =  R  X{u}) 

(4.21) 

P{uj)  =  a^(a;)a(a;) 

(4.22) 

7(0;)  =  a^(a;)X(a;) 

(4.23) 

p{u})  =  X^(a;)X(a>). 

(4.24) 

Then  we  can  rewrite  (4.14)  and  (4.15)  as: 

-  .  .  [R  [(1  -  PM)a(t^)  +  7(w)X(a;)] 

^(a))(l-p(a;))  +  |7(u;)12 

(4.25) 

and 

(  x  Ifi^) 

"^‘^-’-^(a;)(l-p(u;))  +  l7(a;)r 

(4.26) 

whose  implementation  requires  only  the  Cholesky  factorization  of  the  matrix  R  that  is 
independent  of  ui. 


When  applied  to  full  data  sequences,  APES  achieves  the  best  performance  for  a  large 
range  of  filter  lengths  M.  It  can  be  shown  that  increasing  the  filter  length  increases  the 
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resolution  at  the  cost  of  reducing  the  statistical  stability  [7].  Note  that  M  —  N/2  is  the 
maximum  possible  value  of  M  as  beyond  it  R  becomes  singular. 

The  2-D  spectrum  estimation  by  means  of  an  adaptive  filterbank  approach  is  briefly 
explained  as  follows.  Let  denote  an  M  x  M  2-D  impulse  response,  where  M  and  M 

are  the  adaptive  filter  lengths  in  range  and  cross-range,  respectively,  and  let 

h{a;,  (j)  =  vec  [H(w,  w)]  (4-27) 

where  vec[-]  denotes  the  operation  consisting  of  stacking  the  columns  of  a  matrix  on  top 
of  each  other.  Furthermore,  let  X  denote  the  2-D  data  matrix  whose  (n,n)th  element  is 
x{n,  n),  and  Zi  j  denote  the  M  x  M  submatrix  of  X  consisting  of  the  elements  {x{n,  n),n  = 
Z, . . . ,  Z  +  M  -  1,  n  =  r, . . . ,  F-f  M  -  1}  for  Z  =  0, 1, . . . ,  L  -  1  and  F  =  0, 1, . . . ,  L  -  1  where 
L  =  N  —  M  +  1  and  L  =  N  -  M  +  1.  Define 


Zj  r  —  vec[Z;_{], 


(4.28) 


and 


Let 


(4.29) 


1=0  1=0 


a(a),  Lu)  =  a(a))  (S>  a(a;),  (4.30) 

where  <S>  denotes  the  Kronecker  matrix  product,  a(a;)  is  defined  in  (4,8),  and 

a(£D)  =  [l  ...  (4.31) 

Then  the  2-D  adaptive  FIR  filter  and  the  APES  spectral  estimate  can  be  obtained  as  [7]: 


h(u;,  w)  = 


Q”^(a;,a;)a(a;,a;) 


and 


where 


a^(a;,a;)Q”^(a;,a;)a(a;,a;)  ’ 


Q{u,Q)  =  R  -  Z(aj,a))Z^(ci;,a;), 


(4.32) 


(4.33) 


(4.34) 
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with  R  being  the  MM  x  MM  sample  covariance  matrix  associated  with  A  simple 

expression  for  h(a;,a))  can  also  be  obtained  similar  to  (4.25)  by  making  use  of  the  matrix 
inversion  lemma  (cf.  (4.19)). 

Since  the  2-D  APES  algorithm  requires  the  inversion  of  {MM  x  MM)  matrix  (see,  (4.32) 
and  (4.34)),  it  is  computationally  prohibitive  to  apply  the  algorithm  directly  to  data  matrices 
with  large  dimensions.  Instead  we  use  the  chip-based  (approximative)  implementation  of  the 
2-D  APES  algorithm  proposed  in  [25],  as  shown  in  Figure  4.2.  We  apply  the  2-D  EFT  to  X 
to  obtain  an  EFT  image  and  then  break  it  into  small  overlapping  EFT  image  chips  of  size 
Ns  X  Ns  with  Ns  <  N  and  Ns<  N  (in  the  numerical  example,  we  choose  the  overlapping  to 
be  50%  and  set  Ns  =  Ns^  16  while  N  =  N  =  256).  We  perform  the  2-D  IFFT  on  the  small 
chips  to  obtain  the  phase  history  data  chips  to  which  we  apply  the  2-D  APES  algorithm. 
When  applying  the  2-D  APES  algorithm  to  each  phase  history  data  chip,  we  choose  the 
filter  lengths  M  =  Ns/2  and  M  =  iVs/2  and  evaluate  the  spectrum  at  a  frequency  grid  that 
is  4  times  finer  than  the  discrete  FT  (DFT)  grid  in  both  range  and  cross-range.  When  the 
image  chips  are  50%  overlapped,  we  calculate  the  APES  image  chips  Q:(a;,  Q)  only  over  the 
frequency  range  27r^  <  w,  u)  <  27r|  from  the  each  phase  history  data  chip.  The  final  2-D 
APES  SAR  image  consists  of  those  APES  image  chips.  Due  to  overlapping  between  adjacent 
image  chips,  this  procedure  avoids  the  mosaicking  or  tiling  effect,  which  is  noticeable  for  the 
images  shown  in  [26]. 

4.4  GAPES  for  1-D  Gapped  Data 

In  this  section,  we  present  the  algorithm  for  the  spectral  analysis  of  1-D  gapped  data  via 
GAPES.  The  presentation  follows  [27].  The  application  of  the  GAPES  approach  to  the  2-D 
SAR  imaging  with  angle  diversity  data  fusion  will  be  given  in  the  next  section.  We  start  by 
describing  the  algorithm  for  the  spectral  estimation  of  incomplete  data  with  one  gap,  and 
then  extend  it  to  the  data  sequence  with  gaps  of  various  sizes. 
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4.4.1  A  Data  Sequence  with  One  Gap 


In  the  case  that  there  is  only  one  gap  in  the  measured  data  sequence,  we  can  write 


(4.35) 


In  this  case,  we  have  Xa  =  [x^  Xu  =  X2,  and  N  =  Ni  +  N2  +  N3.  The  interpolated  data 
vector  will  be  denoted  by  x  in  which  an  estimate  Xu  of  x^  is  obtained  and  used  to  replace 

Xu. 

We  obtain  initial  APES  estimates  of  h(a;)  and  Q;(a;)  from  the  available  data  Xa  in  the 
following  way.  We  choose  an  initial  filter  length  Mq  such  that  an  initial  full-rank  covariance 
matrix  R  can  be  built  with  the  filter  length  Mq  and  using  Xa  only. 

Let  Lk  —  Nk-Mo  +  1,  k  =  l,3.  We  apply  APES  to  Xa  with  the  following  re-definitions: 

1  /Li— 1  Ni+ATaH-X/s*"! 

^1+^3  \  1=0 


tXi-lU.  .  X 

A(“)  =  E  S(0x''(i)  +  E  x(()x"(l)  (4.37) 

L>1  +  L'S  y  i-Q  1=Ni+N2  ) 

The  APES  spectral  estimate  a:(uj)  as  defined  via  (4.15)  is  a  continuous  function  of  w  and 
hence,  in  principle,  we  could  evaluate  it  at  any  frequency  value.  For  the  initial  estimates 
we  will  evaluate  Qr(tu)  and  h(w)  at  the  DFT  grid:  iVk  —  for  k  =  0, . . . ,  A  1  with 

K  =  N.  (We  can  also  choose  K  >  N,  such  as  K  =  4iV,  for  somewhat  better  results  at  a 
cost  of  more  computations).  The  final  spectral  estimate  can,  if  desired,  be  evaluated  on  a 
finer  grid. 

Next  we  consider  the  estimation  of  Xu  based  on  the  initial  spectral  estimates  &{u))  and 
h(a;)  obtained  as  outlined  above.  Under  the  fairly  natural  assumption  that  the  missing  data 
have  the  same  spectral  content  as  the  available  data  we  can  determine  Xu  from  the  condition 
that  the  output  of  the  filter  h(wfc)  fed  with  the  data  sequence  made  from  Xa  and  Xu  is  as 
close  as  possible  (in  the  LS  sense)  to  a{uk)e^^'‘\  for  i  =  0, . . . ,  L  —  1  and  fc  =  0, . . . ,  A  —  1. 
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Mathematically,  we  propose  to  obtain  Xu  as  the  solution  to  the  following  LS  problem: 

min  Y,  E  |h^(^Jt)xg(0  -  Q;(wfc)e^‘^'='  (4.38) 

k=0  1=0 

where  Xg(0  is  the  overlapping  vector  constructed  from  Xg  in  the  same  way  as  x(l)  is  formed 
from  X  in  (4.6).  Besides  the  intuitive  appeal  of  estimating  x^  in  this  way,  an  additional 
bonus  is  that  we  remain  in  the  APES  framework  (compare  the  criteria  in  (4.13)  and  (4.38)), 
an  observation  to  which  we  will  return  for  further  comments  later  on  in  this  section. 

Let  the  L  x  N  matrix  Hjt  be  defined  by 

h^(wfc) 

Ht=  =1^  ^  3,1  (4.39) 

Ly.Ni  LyN2  LxNs 

h^(Wfc)  _ 

and 

zik  =  d(u;fc)[l  ...  6  (4.40) 

Let  Dfe  =  [Ak  Cfc],  and  define 

y*:  =  Zjt  -  DjtXa.  (4.41) 

Using  this  notation,  we  can  write  the  criterion  in  (4.38)  as: 

min  Y  llBfcXu  -  yikll^  (4-42) 

fc=0 

whose  minimizer  with  respect  to  x^  is  readily  obtained  to  be 

(4.43) 

.fc=0  J  k=0 

Finally,  the  interpolated  data  sequence  x  is  obtained  by  replacing  Xu  with  x^  in  (4.35). 

Once  an  estimate  Xu  has  become  available,  the  next  logical  step  should  consist  of  re- 
estimating  the  spectrum  by  applying  APES  to  the  interpolated  data  sequence  x  made 
from  Xa  and  Xu.  According  to  the  discussion  around  (4.13),  this  entails  the  minimization 
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l2 


(4.44) 


with  respect  to  h(a;it)  and  0!(a;jt)  of  the  function: 


12  |h^K)^(0  - 


it^O  1=0 


subject  to  h^(a;A;)a(cjfc)  =  1,  where  x(Z)  is  made  from  x.  Evidently  the  minimization  of 
(4.44)  with  respect  to  {h(wfc),Q;(a;fc)}f^^  decouples  in  K  minimization  problems  of  the  form 
of  (4.13),  yet  we  prefer  to  write  the  criterion  function  as  in  (4.44)  to  make  the  connection 
with  (4.38).  In  effect,  comparing  (4.38)  and  (4.44)  it  becomes  clear  that  the  alternating 
estimation  of  {a;(a;fc)}  (along  with  {h(a)*:)})  and  Xu  outlined  above  is  nothing  but  a  cyclic 
algorithm  for  solving  the  following  minimization  problem: 


min 

xu,{a(t*;fc),h(a;fc)} 


fc=0  /=0 


(4.45) 


subject  to  h^(ct;fc)a(a;it)  =  1.  As  is  well-known  a  cyclic  minimizer  decreases  the  criterion 
at  each  iteration.  Furthermore,  this  decrease  must  be  strict  in  the  present  case,  as  the 
estimates  of  {Q:(a)fc),h(a;jfe)}  and  Xu  that  we  compute  at  each  iteration  are  uniquely  defined 
(under  very  weak  conditions).  Combining  these  facts  with  the  simple  observation  that  the 
criterion  in  (4.45)  is  nonnegative  and  hence  bounded  from  below,  leads  to  the  conclusion 
that  the  cyclic  minimizer  under  discussion  will  always  converge  to  a  minimum  point  of  the 
criterion.  Whether  that  minimum  point  is  a  local  or  global  one,  this  is  a  harder  question  the 
answer  of  which  will  in  general  depend  on  the  initial  estimate  used. 

The  formulation  in  Equation  (4.45)  has  a  strong  intuitive  appeal  as  it  leads  to: 

(a)  an  analysis  filterbank  {h(a;fc)}  for  which  the  filtered  sequence  is  as  close  as  possible  (in 
the  LS  sense)  to  the  sinusoidal  component  in  x{n)  with  frequency  oJk  ^  [0, 27r), 

(b)  an  estimated  complex  spectrum  that  should  be  almost  leakage-free  (in  view  of  (a))  and 
have  good  statistical  properties  (inherited  from  the  LS  method  used  to  obtain  it),  and 
simultaneously, 

(c)  an  estimate  of  x^  whose  spectral  content  mimics  the  spectral  content  of  the  available 
data  as  much  as  possible  (once  again,  in  the  LS  sense). 

If  we  estimated  Xu  in  a  way  that  were  different  from  (4.38)  then  the  above  cyclic  mini¬ 
mization  interpretation  would  be  lost  (this  interpretation  is  the  bonus  of  estimating  x^  in  the 
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APES  framework,  which  was  mentioned  following  (4.38)).  The  same  thing  would  happen  if 
we  estimated  the  spectrum  of  the  interpolated  sequence  by  a  method  that  is  different  from 
APES. 

The  GAPES  algorithm  implements  the  cyclic  minimization  of  (4.45)  outlined  above. 

A  step-by-step  summary  of  this  algorithm  follows. 

Step  0:  Initialization:  start  the  algorithm  using  the  initialization  approach  mentioned  at 
the  beginning  of  this  section. 

Step  1:  Gap  filling:  update  the  gap  estimate  in  {x{n)}  with  Xu  obtained  by  using  (4.43) 
based  on  the  most  recently  determined  adaptive  FIR  filter  and  APES  spectrum  {h(a;fc),  a(a;fc)}. 
Step  2:  Spectrum  estimation:  determine  the  adaptive  FIR  filter  and  APES  spectrum  with 
(4.25)  and  (4.26),  respectively,  from  the  most  recently  interpolated  data  sequence  {f  (n)}. 
Step  3:  Iteration:  repeat  Steps  1  and  2  till  the  relative  change  of  the  cost  function  in  (4.45) 
is  equal  to  or  less  than  a  preset  threshold  value  e  or  the  number  of  iterations  reaches  a  preset 

number  /max- 

4.4.2  Spectral  Estimation  of  1-D  Data  Sequences  with  Multiple  Gaps 

The  extension  of  the  algorithm  presented  in  the  previous  section  to  the  multi-gap  case 
is  straightforward.  One  approach  could  be  to  let  the  gap  parameter  vector  Xu  contain  all 
gaps,  accordingly  modify  (4.42)  and  all  related  equations,  and  then  obtain  an  estimate  of 
missing  data  simultaneously.  However,  in  many  applications,  the  measured  data  can  be 
non-stationary  but  locally  stationary.  Therefore,  we  intend  to  perform  the  gap  filling  for 
each  individual  gap  separately.  Each  gap  is  filled  in  by  making  use  of  only  the  two  adjacent 
continuous  segments. 

4.5  SAR  Imaging  for  Angle  Diversity  Data  via  GAPES 

The  1-D  GAPES  algorithm  can  be  readily  applied  to  SAR  imaging  with  angle  diversity 
data  fusion,  where  the  2-D  phase  history  data  matrix  Xg  whose  (n,  n)th  element  is  Xg(n,  n) 
has  missing  columns.  The  principal  steps  are  summarized  as  follows. 
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Step  I:  Perform  1-D  WFFT  for  each  available  column  of  the  phase  history  data  matrix  Xg 
to  obtain  Yi. 

Step  II:  Fill  in  the  gaps  for  each  row  of  Yi  via  GAPES  to  obtain  Y2. 

Step  III:  If  a  2-D  WFFT  SAR  image  is  desired,  then  compute  1-D  WFFTs  for  each  row  of 
Y2  and  stop.  Otherwise,  go  to  the  next  step. 

Step  IV:  Apply  1-D  IFFT  and  de-windowing  to  each  column  of  Y2  to  obtain  the  gap-filled 
phase  history  data  matrix  X. 

Step  V:  Obtain  the  2-D  SAR  image  by  applying  the  chip-based  2-D  APES  outlined  above 
to  the  data  matrix  X. 

4.6  Numerical  Results 

We  now  present  some  numerical  examples  to  illustrate  the  performance  of  the  GAPES 
algorithm  for  spectral  analysis  of  gapped  data  and  SAR  imaging  with  angle  diversity  data 
fusion.  We  choose  K  —  N  ior  the  iteration  steps  and  the  spectrum  is  estimated  at  discrete 
frequencies  a;*  =  2'Kk/K,  or  fk  =  Wfe/27r,A:  =  O,---,  AT  -  1.  The  simulated  signals  in  the 
examples  are  corrupted  by  additive  zero-mean  white  Gaussian  noise  with  variance  0.01.  In 
the  initialization  step  for  each  gap,  the  filter  length  is  chosen  as  the  smallest  integer  larger 
or  equal  to  3/4  of  the  length  of  the  shorter  adjacent  continuous  segment.  We  set  c  =  10~^ 
to  detect  the  convergence  of  the  algorithm  and  the  maximum  number  of  iteration  is  set  to 
Imax  =  10-  After  the  estimation  procedure  terminates,  we  obtain  a  finer  spectral  estimate 
with  APES  by  using  K  ■=  AN  for  the  1-D  spectral  estimation  example. 

We  compare  our  algorithm  with  the  WFFT  and  with  the  parametric  CLEAN  algorithm 
(see  e.g.  [16]).  1-D  and  2-D  Taylor  windows  with  order  5  and  sidelobe  level  -35  dB  are  used 
for  1-D  and  2-D  WFFTs,  respectively.  When  using  CLEAN,  we  model  each  scatterer  of  the 
target  of  interest  as  a  complex  sinusoid  with  constant  amplitude  and  phase.  CLEAN  first 
estimates  the  parameters  of  the  strongest  scatterer,  subtracts  it  from  the  original  signal,  and 
then  repeats  this  process  for  the  next  strongest  scatterer  until  the  predetermined  number 
of  scatterers  is  reached  or  the  estimated  amplitude  of  the  current  scatterer  is  small  enough 
compared  to  the  strongest  one.  Once  the  scatterer  parameters  are  estimated,  they  are  used 
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to  simulated  the  data  in  the  gaps.  The  WFFT  approach  is  then  used  to  form  the  CLEAN 
spectrum. 

We  also  consider  the  Autoregressive  (AR)  Model  based  Extrapolation  approach,  referred 
to  as  ARME.  Each  available  data  segment  (except  for  the  last  one)  is  used  with  the  Yule- 
Walker  method  [22,  23,  24],  which  guarantees  stability  for  forward  linear  prediction  to  fill  in 
one  half  of  the  gap  to  the  right.  Each  available  data  segment  (except  for  the  first  one)  is  also 
used  with  the  Yule- Walker  method  for  backward  linear  prediction  to  fill  in  one  half  the  gap 
to  the  left.  We  choose  the  AR  model  order  to  be  one  half  of  the  length  of  the  available  data 
segment  when  it  is  used  to  fill  in  one  half  of  the  gap  to  the  right  or  to  the  left  of  the  data 
segment.  The  WFFT  approach  is  then  used  to  form  the  ARME  spectrum  after  all  the  gaps 
are  filled.  We  remark  that  although  the  Yule- Walker  method  avoids  the  instability  problem 
for  data  extrapolation,  the  extrapolated  data  sequences  exponentially  delay  to  zero,  which 
can  cause  artifacts  (see  the  examples  below). 

Figure  4.3  presents  an  example  of  spectral  analysis  of  incomplete  data  with  two  gaps  of 
different  sizes.  The  total  length  of  the  data  sequence  is  AT  =  128  and  samples  23  through  46 
and  76  through  100  are  missing  (note  that  almost  40  percent  of  the  data  is  missing).  Figure 
4.3(a)  shows  the  true  spectrum  which  consists  of  four  spectral  lines  located  at  fi  —  0.05, 
/2  =  0.065,  fi  =  0.26,  and  /z  =  0.28  with  complex  amplitudes  ai  =  012  =  0:3  =  1  and 
a4  =  0.5,  respectively,  and  a  continuous  spectral  component  centered  at  /«  =  0.18  with  a 
width  b  =  0.015  and  a  constant  modulus  a  =  0.25.  In  Figure  4.3(b)  the  WFFT  is  applied 
to  the  data  by  filling  in  the  gaps  with  zeros  (note  the  artifacts  in  the  spectrum  due  to  the 
missing  data).  Figures  4.4(c)  and  (d)  show  the  moduli  of  the  WFFT  and  APES  spectra  of 
the  complete  data  sequence  (note  the  superior  resolution  of  APES  as  compared  to  WFFT). 
Figures  4.3(e)  and  (f)  demonstrate  the  moduli  of  the  WFFT  spectra  of  the  interpolated 
data  sequences  via  CLEAN  and  ARME,  and  Figures  4.3(g)  and  (h)  present  the  moduli  of 
the  WFFT  and  the  APES  spectra  of  the  data  sequence  interpolated  via  GAPES.  Note  that 
in  this  case,  CLEAN  and  ARME  cannot  be  used  to  effectively  eliminate  the  spectral  artifacts 
due  to  the  amount  of  data  missing  whereas  GAPES  is  shown  to  be  effective  for  filling  in  the 
gaps  and  estimating  the  spectrum.  In  Figure  4.3(i),  we  show  the  true  data  and  the  data 
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interpolated  via  GAPES  (only  the  real  part  is  shown;  the  behavior  of  the  imaginary  part  is 
similar).  The  results  displayed  so  far  are  for  one  randomly  picked  realization  of  the  data. 
In  Figure  4.3(j),  we  show  the  GAPES  spectrum  for  five  different  randomly  selected  data 
realizations. 

We  now  illustrate  the  performance  of  SAR  imaging  with  angle  diversity  data  fusion  via 
GAPES  with  the  high  resolution  phase  history  data  of  a  Slicy  object  at  0°  azimuth  angle 
generated  by  XPATCH  [28],  a  high  frequency  electromagnetic  scattering  prediction  code  for 
complex  3-D  objects.  A  photo  of  the  Slicy  object  taken  at  45°  azimuth  angle  is  shown  in 
Figure  4.4(a).  The  original  XPATCH  data  matrix  has  a  resolution  of  0.043  meters  in  both 
range  and  cross-range  and  a  size  oi  N  —  N  —  256.  Figures  4.4(b)  and  (c)  show  the  2-D 
WEFT  and  2-D  chip-based  APES  SAR  images  from  the  complete  XPATCH  data  matrix, 
respectively. 

To  demonstrate  the  advantage  of  the  angle  diversity  data  fusion  via  GAPES,  we  artifi¬ 
cially  create  gaps  as  follows.  We  consider  the  first  and  the  last  8  columns  of  the  XPATCH 
data  matrix  as  missing  and  let  the  remaining  240  columns  contain  7  uniformly  spaced  gaps 
each  with  16  columns.  Each  continuous  data  segment  contains  16  columns  as  well.  The 
overall  missing  data  percentage  is  50%.  Figure  4.5(a)  exhibits  the  WEFT  image  by  filling  in 
the  gaps  with  zeros.  Note  that  in  cross-range,  the  point-like  scatterers  are  split  into  multiple 
separate  point  scatterers  due  to  the  missing  data  so  that  distinct  ghosts  appear.  Figures 
4.5(b)  and  (c)  show  the  SAR  images  obtained  by  using  2-D  WEFT  on  the  CLEAN  and 
ARME  interpolated  data,  respectively.  Shown  in  Figure  4.5(d)  is  the  SAR  image  obtained 
by  using  2-D  WFFT  on  the  GAPES  interpolated  data  while  the  final  GAPES  image  ob¬ 
tained  by  using  the  2-D  chip-based  APES  on  the  GAPES  interpolated  data  is  presented  in 
Figure  4.5(e).  Note  that  the  spectral  lines  due  to  the  dihedrals  are  distorted  by  CLEAN  due 
to  the  inherent  modeling  of  them  as  point  scatterers.  The  ARME  approach  cannot  over¬ 
come  the  missing  data-induced  ghost  problem  associated  with  the  point-like  scatterer  due 
to  the  linear  prediction  model  not  matching  the  sinusoidal  signal  for  the  point-like  scatterer. 
From  Figures  4.4  and  4.5,  we  observe  that  the  SAR  images  obtained  by  using  the  GAPES 
interpolated  data  are  almost  the  same  as  those  obtained  by  using  the  complete  data. 
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4.7  Conclusions 


We  have  investigated  the  application  of  the  APES  (Amplitude  and  Phase  Estimation) 
filter  to  spectral  analysis  of  1-D  gapped  data  sequences  and  for  2-D  SAR  imaging  with  angle 
diversity  data  fusion.  The  proposed  GAPES  (Gapped-data  APES)  algorithm  iterates  the 
following  two  steps:  (1)  estimating  the  adaptive  filter  and  the  corresponding  spectrum  via 
APES,  and  (2)  filling  in  the  gaps  via  a  least  squares  (LS)  fitting.  For  the  SAR  imaging 
with  angle  diversity  data  fusion,  we  performed  the  1-D  WFFTs  in  range,  used  the  GAPES 
algorithm  to  interpolate  the  gaps  in  the  aperture  for  each  range,  applied  the  1-D  IFFTs  and 
de-windowed  in  range,  and  applied  the  2-D  APES  algorithm  to  the  interpolated  matrix  to 
obtain  the  final  2-D  SAR  image.  Numerical  examples  have  been  presented  to  demonstrate 
the  effectiveness  of  the  algorithm  and  compare  it  with  the  windowed  FFT  approach,  the 
parametric  CLEAN  method,  and  the  ARME  (AR  Model  based  Extrapolation)  approach. 
The  GAPES  algorithm  has  been  shown  to  provide  the  best  spectral  estimates  and  SAR 
images. 
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Figure  4.1:  Schematic  map  of  a  spotlight  SAR  data  acquisition  with  angle  diversity. 


Figure  4.2:  Breaking  a  large  data  matrix  into  small  chips  and  applying  2-D  APES  to  each 
small  chip  to  obtain  the  2-D  spectral  estimation  of  a  data  matrix  with  large  dimensions. 


Figure  4.3:  Moduli  of  complex  spectral  estimates  (solid  lines  in  (b)-(h)  compared  with  the 
true  spectrum  (dotted  lines  in  (b)-(h));  (a)  true  spectrum,  (b)  WFFT  spectrum  of  the 
gapped  data  (filling  in  the  gaps  with  zeros),  (c)  WFFT  spectrum  of  the  complete  data,  (d) 
APES  spectrum  of  the  complete  data,  (e)  WFFT  spectrum  of  the  interpolated  data  sequence 
via  CLEAN,  (f)  WFFT  spectrum  of  the  interpolated  data  sequence  via  ARME,  (g)  WFFT 
spectrum  of  the  interpolated  data  sequence  via  GAPES,  and  (h)  APES  spectrum  of  the 
interpolated  data  sequence  via  GAPES,  (i)  True  data  and  data  interpolated  via  GAPES, 
and  (j)  GAPES  spectrum  for  five  different  realizations. 
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(a) 


Figure  4.4;  (a)  Target  photo  taken  at  45“  azimuth  angle,  and  (b)  2-D  WFFT  (c)  and  2-D 
APES  SAR  images  obtained  by  using  the  complete  XPATCH  data  of  size  256x256. 
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Figure  4.5:  Comparison  of  SAR  images,  (a)  WFFT  image  of  gapped  data  matrix,  (b)  WFFT 
image  with  angle  diversity  data  fusion  via  CLEAN,  (c)  WFFT  image  with  angle  diversity 
data  fusion  via  ARME,  (d)  WFFT  image  with  angle  diversity  data  fusion  via  GAPES,  and 
(e)  APES  image  with  angle  diversity  data  fusion  via  GAPES. 
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5.  A  Quasi-parametric  Algorithm  for  SAR  Target 
Feature  Extraction  and  Imaging  with  Angle  Diversity 


5.1  Introduction 

Target  feature  extraction  from  spotlight  mode  synthetic  aperture  radar  (SAR)  measure¬ 
ments  plays  an  important  role  in  many  applications  including  battlefield  awareness  [1]  and 
automatic  target  recognition  (ATR)  [2,  3].  Fast  Fourier  transform  (FFT)  based  approaches 
have  been  widely  used  due  to  their  high  computational  efficiency  and  robust  performance. 
However,  the  disadvantages  of  such  approaches  are  their  low  resolution,  poor  accuracy,  and 
high  sidelobes.  A  variety  of  spectral  estimation  methods  have  been  proposed  to  obtain  target 
features  with  high  resolution  and  low  sidelobes.  Of  these,  the  parametric  (e.g.  [4,  5,  6,  7]), 
non-parametric  (e.g.  [8,  9,  10,  11]),  and  semi-parametric  [12]  approaches  have  been  con¬ 
sidered  for  SAR  image  formation  and  target  feature  extraction.  In  general,  parametric 
approaches  may  outperform  their  nonparametric  counterparts  in  resolution  and  accuracy 
but  are  more  sensitive  to  data  modeling  errors. 

Azimuth  angle  diversity  data  fusion,  which  combines  several  sets  of  radar  measurements 
made  at  different  azimuth  angle  locations  with  advanced  signal  processing  methodology, 
provides  a  useful  tool  for  SAR  target  feature  extraction  and  imaging  with  improved  resolution 
and  can  be  used  to  achieve  better  ATR  performance.  Prom  the  viewpoint  of  signal  processing, 
the  solution  to  the  angle  diversity  data  fusion  problem  shares  a  similar  principle  to  the  missing 
or  gapped  data  problem,  which  occurs,  for  example,  in  astronomy.  Valuable  astronomical 
data  are  often  available  only  in  the  form  of  groups  of  samples  which  are  gapped  within 
rather  long  intervals  when  no  reliable  measurements  can  be  taken  (see,  e.g.  [13,  14,  15,  16], 
and  the  references  therein).  Special  attentions  have  been  paid  to  detecting  the  presence  of 
one  or  more  periodic  signals  from  incomplete  measurements  using  techniques  such  as  the 
periodogram  approach  [14],  the  CLEAN  deconvolution  [15],  and/or  reconstruction  of  the 
unevenly  sampled  data  on  a  regular  grid  by  making  certain  assumptions  on  the  data. 

Most  parametric  SAR  target  feature  extraction  methodologies  make  the  assumption  that 
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the  man-made  target  of  interest  consists  of  several  trihedral  corner  reflectors  (point-like 
scatterers)  and  derive  the  corresponding  algorithms  based  on  the  two-dimensional  (2-D) 
complex  sinusoidal  data  model  in  which  the  amplitude  and  phase  are  constant  in  both  range 
and  cross-range.  However,  this  data  model  does  not  always  flt  the  actual  radar  observations. 
Many  man-made  targets  can  contain  dominant  scattering  mechanisms  other  than  trihedral- 
like  scatterers.  For  instance,  radar  responses  from  many  man-made  targets  such  as  vehicles 
and  buildings  are  primarily  caused  by  both  trihedral  and  dihedral  corner  reflectors  [17].  A 
principal  difference  between  a  trihedral  and  a  dihedral  corner  reflector  lies  in  their  radar 
responses  in  cross-range.  That  is,  the  former  can  be  modeled  as  a  complex  sinusoid  with  a 
constant  amplitude  and  phase,  while  the  latter  can  be  approximately  described  as  a  complex 
sinusoid  with  its  amplitude  being  a  sine  function  (deflned  as  sinc(a:)=sin(a:)/a;)  and  its  phase 
a  constant.  A  mixed  data  model  was  presented  in  [7],  which  characterizes  the  reflections 
in  cross-range  from  a  trihedral  and  a  dihedral  corner  reflector  with  a  constant  and  a  sine 
function,  respectively.  In  [12],  this  data  model  is  extended  by  modeling  each  target  scatterer 
as  a  2-D  complex  sinusoid  with  arbitrary  amplitude  and  constant  phase  in  cross-range  and 
with  constant  amplitude  and  phase  in  range,  and  a  Semi-PARametric  (SPAR)  algorithm 
was  proposed  for  SAR  target  feature  extraction  and  high  resolution  image  formation.  Due 
to  its  flexible  data  model,  SPAR  combines  the  advantages  of  both  parametric  and  non- 
parametric  spectral  estimation  methods,  and  thus  has  batter  estimation  performance  and 
higher  resolution  capability  than  nonparametric  algorithms  and  is  more  robust  against  data 
modeling  errors  over  parametric  ones. 

To  achieve  excellent  ATR  performance,  the  radar  range  resolution  capability  is  becoming 
an  increasingly  important  factor  to  be  considered  in  the  radar  system  design.  For  example, 
ultra  wideband  (UWB)  radar  has  been  an  effective  technology  due  to  its  extremely  large 
bandwidth  (several  GHz)  and  very  high  range  resolution  (on  the  order  of  several  centimeters). 
As  the  range  resolution  increases,  the  radar  cross-sections  of  target  scatterers  vary  within 
the  large  bandwidth  [18],  and  the  assumption  of  constant  responses  can  no  longer  hold. 
This  makes  the  signal  processing  needed  by  such  ultra  high  range  resolution  radars  more 
difficult.  More  flexible  data  model  is  required  to  achieve  robust  performance  of  target  feature 
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extraction. 

In  this  paper,  we  study  the  SAR  target  feature  extraction  and  imaging  with  angle  diversity 
using  an  UWB  radar.  Instead  of  having  a  radar  collecting  data  continuously  over  a  large  angle 
to  achieve  high  resolution  in  cross-range,  which  is  a  huge  commitment  of  radar  resources,  we 
attempt  to  achieve  high  resolution  in  cross-range  by  collecting  data  continuously  over  a  small 
angle  but  at  several  diverse  azimuth  angle  locations.  More  specifically,  with  angle  diversity, 
the  radar  collects  data  continuously,  stops,  collects  more  data  continuously,  stops  again,  and 
so  forth.  Hence  the  data  available  with  angle  diversity  consist  of  several  pieces  of  the  data 
collected  continuously  over  a  large  angle.  Instead  of  using  the  mixed  data  model  in  [7]  or  the 
semi-parametric  data  model  in  [12],  we  use  a  more  flexible  data  model,  which  describes  each 
target  scatterer  as  a  2-D  complex  sequence  with  arbitrary  amplitude  and  constant  phase 
in  both  range  and  cross-range.  This  data  model  is  essentially  quasi-parametric  because  of 
the  arbitrary  amplitude  assumption  in  both  dimensions.  The  feature  extraction  algorithm 
based  on  such  a  flexible  data  model  is  more  robust  against  the  data  modeling  errors  than 
the  parametric  and  semi-parametric  methods.  A  new  algorithm,  referred  to  as  the  QUALE 
(QUasi-parametric  ALgorithm  for  target  feature  Extraction)  algorithm,  is  presented  for  the 
SAR  target  feature  extraction  with  angle  diversity.  QUALE  first  chips  each  scatterer  out 
from  the  image  of  the  target  of  interest  and  apply  2-D  inverse  FFT  (IFFT)  to  obtain  the 
phase  history  data.  Secondly  QUALE  estimates  the  model  parameters  involved  in  the  flexible 
data  model,  which  includes,  for  each  scatterer  of  the  target  of  interest,  locations  in  range 
and  cross-range,  a  constant  phase,  and  a  2-D  arbitrary  real-valued  amplitude  sequence. 
Thirdly  QUALE  performs  an  average  in  the  range  dimension  over  the  estimated  2-D  real- 
valued  amplitude  sequence,  models  the  so-obtained  1-D  sequence  with  a  simple  sine  function, 
and  estimates  the  relevant  sine  parameters  by  minimizing  a  nonlinear  least  squares  (NLS) 
function.  Finally,  the  2-D  SAR  image  is  reconstructed  by  using  the  estimated  features.  Like 
the  semi-parametric  data  model  in  [12],  there  are  also  ambiguity  problems  associated  with 
our  data  model.  To  alleviate  this  problem,  QUALE  uses  an  isolation  process  to  chip  out  the 
most  dominant  scatterer  with  a  2-D  window  in  the  image  domain,  so  that  it  only  deals  with 
one  scatterer  at  a  time. 
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The  remainder  of  this  paper  is  organized  as  follows.  Section  2  establishes  the  data  model 
and  formulates  the  problem  of  interest.  The  QUALE  algorithm  is  presented  in  Section  3.  In 
Section  4,  numerical  examples  are  presented  to  illustrate  the  performance  of  the  proposed 
algorithm.  Finally,  Section  5  contains  our  conclusions. 

5.2  Data  Model 

To  obtain  high  resolution  SAR  target  features  and  to  reconstruct  the  SAR  images  of 
targets  of  interest,  we  need  to  build  a  proper  data  model  for  the  target  scatterers.  It  has 
been  investigated  in  [12]  that  the  assumption  of  arbitrary  amplitude  and  constant  phase 
in  cross-range  is  in  general  more  applicable  than  the  mixed  data  model  to  characterize  the 
scattering  mechanisms  including  both  trihedral  and  dihedral  corner  reflectors.  Since  for  an 
UWB  radar  collecting  data  continuously  over  a  large  angle,  the  RCS  of  a  target  scatterer 
is  not  constant  within  the  radar  bandwidth,  we  model  the  received  signal  reflected  from  a 
target  scatterer  as: 

s(n,n)  =  n  =  0,l,---,Ar- 1,  n  =  0, 1,  •  •  • ,  AT  -  1,  (5.1) 

where  N  and  N  denote  the  dimensions  of  data  samples  in  range  and  cross-range,  respectively; 
z{n,  n)  is  an  arbitrary  unknown  2-D  real- valued  function  of  n  and  n  determined  by  the  RCS 
of  the  scatterer  as  a  function  of  frequency  and  angle  with  respect  to  the  radar  and  the 
particular  scattering  mechanism  of  the  scatterer;  ^  is  a  constant  phase;  finally,  {/,  /}  is  the 
frequency  pair  proportional  to  the  locations  in  range  and  cross-range  of  the  scatterer.  This 
data  model  is  essentially  quasi-parametric  since  little  parameterization  is  assumed  in  either 
range  or  cross-range. 

Assume  that  a  target  of  interest  consists  of  K  scatterers.  Then  the  target  data  model  in 
the  presence  of  noise  and  clutter  has  the  form: 

K 

y{n,n)  =  +  e(n,h), 

fc=i 

n  =  0,l,---,Ar-l,  n  =  0,l,---,iV-l,  (5.2) 
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where  Zk{n,n)  denotes  an  arbitrary  2-D  real- valued  amplitude  function  of  n  and  n  for  the 
kth  scatterer;  (pk  and  {/*;,  /*:},  respectively,  are  the  constant  phase  and  the  frequency  pair 
of  the  A:th  scatterer;  finally,  e(n,  n)  denotes  the  unknown  2-D  noise  and  clutter  sequence. 
Note  that,  if  the  amplitude  function  Zk{n,n)  is  assumed  to  be  constant  in  range,  we  get  the 
same  data  model  as  the  one  presented  in  [12]. 

For  angle  diversity  data  fusion,  the  data  measurements  are  recorded  intermittently.  For 
this  case,  the  data  model  can  be  expressed  as: 


y(n,  n) 


! 


0, 

y(n,n). 


71  =  0, 1,  •  •  • ,  iV  —  1,  tI  G  G, 
otherwise, 


(5.3) 


where  G  consists  of  the  indices  where  the  data  measurements  are  not  made.  Note  that,  (5.3) 
is  equivalent  to  setting 


and  letting 


Zk{n,  n) 


0, 

5fc(n,n), 


71  =  0,  1,  •  ■  • ,  iV  —  1,  71  G  G, 

otherwise, 


K 


y{n,  n)  =  ^  Zk{n,  7i)e’^e^2ir(/fcn+/tn)  _j_ 

A:=l 

71  =  0, 1,  •  •  • ,  iV  —  1,  fi  =  0, 1,  •  •  * ,  iV  —  1, 


(5.4) 


(5.5) 


where  e{n,  n)  is  formed  from  e{n,  n)  in  the  same  way  as  Zk{n,  n)  is  formed  from  Zk{n,  n).  To 
reduce  the  amount  of  computations  and  avoid  model  ambiguity,  we  chip  out  each  scatterer 
one  at  a  time  (see  the  next  section  for  more  details).  Our  objective  is  to  estimate  the  model 
parameters  from  each  scatterer  chip. 

Once  the  amplitude  estimate  of  each  scatterer  is  obtained,  we  attempt  to  obtain  approx¬ 
imate  target  features  by  averaging  it  over  n.  We  remark  that  we  can  readily  modify  our 
approach  to  model  the  amplitude  estimate  to  vary  with  ti  as  7i“  [18].  We  choose  to  average 
it  over  n  to  simplify  the  problem. 

We  then  approximately  model  the  /cth  real-valued  amplitude  average  as  a  sine  function 
of  n,  i.e.,  as  arsine  [7r6fc(n  -  Tjt)],  where  bk,  and  Tk  denote  the  unknown  real- valued  am¬ 
plitude,  the  spectral  width,  and  the  signal’s  peak  location  of  the  kth  scatterer,  respectively. 
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These  three  parameters  have  quite  clear  physical  interpretations.  That  is,  cuk  is  proportional 
to  the  RCS,  bk  is  related  to  the  length  of  the  scatterer  in  cross-range,  and  Tk  gives  the  peak 
location  of  the  data  sequence  and  is  determined  by  the  orientation  of  the  scatterer.  We 
estimate  the  approximate  scatterer  parameters  {ak,  bk,  r*}  from  the  amplitude  average  and 
reconstruct  an  approximate  2-D  SAR  image  using  the  estimated  parameters. 

5.3  The  QUALE  Algorithm 

The  QUALE  algorithm  consists  of  three  main  steps.  First,  each  scatterer  is  chipped  out 
in  the  image  domain  by  using  an  isolation  method,  which  will  be  described  in  Section  5.3.1, 
then  2-D  IFFT  is  applied  to  the  image  chip  to  get  the  corresponding  phase  history  data  of 
each  scatterer.  Secondly,  the  model  parameters  of  each  scatterer  are  estimated  according  to 
the  2-D  quasi-parametric  model  in  (5.5)  (see  Section  5.3.2  for  details),  then  the  sine  function 
parameters  are  estimated  by  minimizing  a  NLS  cost  function  (see  Section  5.3.2  for  details). 
Finally,  the  phase  history  data  of  each  scatterer  are  simulated  according  to  the  estimated 
scatterer  parameters  and  2-D  FFT  is  applied  to  the  simulated  data  to  obtain  the  enhanced 
image  chips,  and  then  these  enhanced  image  chips  are  put  back  to  the  locations  of  the  original 
image  chips  to  reconstruct  the  2-D  approximate  SAR  image  (see  Section  3.3  for  details). 

5.3.1  Scatterer  Isolation 

Because  of  the  flexibility  of  the  data  model  we  use,  there  will  be  ambiguity  problems  when 
two  or  more  scatterers  are  located  in  the  same  range  or  cross  range.  To  avoid  this  problem, 
we  chip  each  scatterer  out  and  work  only  on  one  of  them  at  a  time.  Before  we  present  the 
isolation  method,  we  need  to  consider  another  problem:  when  working  with  the  gapped  data 
with  angle  diversity,  the  peak  of  a  single  scatterer  will  be  split  into  several  peaks,  which 
makes  it  difficult  to  single  out  each  scatterer.  This  phenomenon  is  illustrated  in  Figures 
5.1(a)  and  (b).  Figure  5.1(a)  is  the  zero-padded  FFT  of  an  all-one  sequence,  and  Figure 
5.1(b)  is  the  zero-padded  FFT  of  the  gapped  data  due  to  angle  diversity.  More  specifically, 
the  length  of  the  all-one  sequence  is  32,  the  gapped  data  due  to  angle  diversity  is  formed  by 
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setting  the  middle  50%  to  zero.  The  two  sequences  are  both  zero-padded  to  128.  It  is  clear 
that  one  peak  is  split  into  many  peaks.  We  avoid  this  peak  splitting  problem  by  isolating  the 
scatterers  based  on  the  image  generated  from  a  continuous  data  segment  that  has  the  most 
power.  Let  X  denote  the  original  data  matrix  and  X  denote  the  continuous  data  segment 
with  the  most  power,  with  both  of  which  N  x  N  matrices  (X  can  be  obtained  by  setting 
other  data  segments  to  zero).  Let  V  denote  the  2-D  FFT  of  X  without  zero-padding,  and 
V  denote  the  2-D  FFT  of  X  with  zero-padding  to  Lx  L,  where  zero-padding  is  necessary  to 
improve  the  accuracy.  We  determine  the  size  and  position  of  the  isolation  window  of  each 
scatterer  on  V  and  chip  out  the  scatterer  from  V. 

There  is  a  linear  relationship  on  the  size  and  position  between  the  windows  in  V  and  V. 
If  {1, 1}  is  a  location  in  V,  then  the  corresponding  location  {n,  n}  in  V  can  be  obtained  with 

"  =  [f  ('  -  1)  +  ij  .  (5-6) 

and, 

fl=  [  j(/-l)  +  lJ,  (5.7) 

where  [-J  denotes  the  nearest  integer  towards  minus  infinity.  Next  we  describe  how  to 
determine  the  windows  needed  to  chip  out  each  scatterer.  In  [12],  a  windowing  method  is 
proposed  to  chip  out  the  most  dominant  scatterer  one  at  a  time.  In  this  paper,  we  use  a  more 
robust  version  to  meet  the  requirement  of  the  gapped  data  due  to  angle  diversity.  Unlike 
SPAR,  which  uses  the  same  threshold  value  for  both  range  and  cross-range,  we  use  different 
threshold  values  for  range  and  cross-range. 

To  chip  out  the  most  dominant  scatterer,  we  first  determine  the  peak  location 
from  the  magnitude  of  V.  We  then  fix  /  to  f''  and  search  around  by  starting  from  to 
determine  the  interval  li  <  I  <  h  so  that  the  magnitude  of  V  is  above  a  certain  threshold 
Tf.  Similarly,  we  can  fix  I  to  l'^  and  search  around  by  starting  from  L*"  to  determine  the 
interval  Ii  <  I  <12  so  that  the  magnitude  of  V  is  above  another  certain  threshold  Ter-  The 
corresponding  window  position  {711,712,  ni,h2}  on  V  can  be  determined  according  to  (5.6) 
and  (5.7).  The  corresponding  phase  history  data  of  the  scatterer.  can  be  obtained  by  chipping 
out  the  image  chip  determined  by  {tii,  712,  m,  712}  from  V  and  applying  2-D  IFFT.  Note  that 
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when  two  scatterers  are  closely  spaced  in  range,  they  may  fall  into  the  same  window,  which 
can  result  in  biased  target  feature  estimates.  To  minimize  the  problem  while  determining 
the  interval  k  <!<  k,  we  also  check  the  ratio  R  between  the  height  of  the  second  peak  and 
that  of  the  strongest  peak  within  the  window  around  {l"^,  f"^).  If  it  is  larger  than  a  threshold 
Tratio,  they  will  be  split  into  two  separate  windows  with  the  valley  point  between  the  two 
peaks  as  the  window  boundary. 

The  steps  of  chipping  out  all  of  the  scatterers  can  be  summarized  as  follows: 

Step  0;  Compute  the  2-D  FFT  V  and  V  of  X  and  X,  respectively. 

Step  k:  Determine  the  window  Wk  for  the  strongest  scatterer  on  V,  compute  the  corre¬ 
sponding  window  Wk  on  V.  Set  those  elements  of  V  covered  by  Wk  to  zero.  Obtain  the 
corresponding  phase  history  data  matrix  for  the  scatterer  by  applying  2-D  IFFT  to  the  image 
chip  obtained  by  applying  Wk  to  V. 

Stop:  Stop  when  we  have  K  scatterers,  which  can  be  predetermined  or  according  to  how 
strong  the  new  scatterer  is  as  compared  to  the  strongest  scatterer. 

5.3.2  Scatterer  Parameter  Estimation 

We  first  present  the  estimation  of  the  parameters  of  each  scatterer  according  to  the  2-D 
quasi-parametric  model,  then  the  sine  function  parameters  are  estimated  by  minimizing  a 
NLS  cost  function. 

2-D  Quasi-Parametric  Parameter  Estimation 

Due  to  the  scatterer  isolation  process  described  in  Section  3.1,  the  phase  history  data 
matrix  corresponding  to  the  A:th  scatterer  in  the  presence  of  noise  has  the  form: 

yk{n,  n)  =  Xk{n,  ek{n,  n), 

n  =  0, 1,- •  •  —  1,  n  =  0,l,-*-,iVfc  -  1,  (5.8) 

where  Nk  and  Nk  are  the  dimensions  of  the  fcth  scatterer  data  matrix  in  range  and  cross¬ 
range,  respectively;  Xk{n,  h)  denotes  an  arbitrary  2-D  real- valued  amplitude  function  of  n  and 
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n  for  the  scatterer;  and  {/*,  /t},  respectively,  are  the  constant  phase  and  the  frequency 
pair  of  the  scatterer;  finally,  ejt(n,  n)  denotes  the  unknown  2-D  noise  and  clutter  sequence. 
Let  Yjfc,  Xfc,  and  be  Nk  x  Nk  matrices  with  their  (n,  n)th  elements  being  yk{n,  n),  Xk{n,  n) 
and,  eit(n,n),  respectively.  Let  T>{fk,fk)  denote  the  following  Nk  x  Nk  matrix; 

^{fk,  Ik)  =  (5.9) 

where  (•)^  denotes  the  transpose, 

=  •••  (5.10) 

and 

<.'R.(/0  =  [l  (5.11) 

Then  (5.8)  can  be  rewritten  in  a  matrix  form  as 

Yfc  =  e^'^*D(/)t,  Ik)  O  Xfc  +  Efc,  (5.12) 

where  ©  denotes  the  Hadamard  matrix  product,  i.e.,  the  element-wise  matrix  product.  The 
NLS  estimates  {Xfc,  ^k,  Ik,  h}  of  {Xfc,  <f>k,  fk,  Ik},  can  be  obtained  by  minimizing  the  follow¬ 
ing  cost  function 

Ci{Xt,h,fkJk)  =  ||Yi-e»*D(A,/k)©Xtl|^ 

Nk-l  _  2 

n=0  n=:0 

where  ||  •  Hi?  denotes  the  Frobenius  norm.  After  some  straightforward  derivations,  we  can 
rewrite  (5.13)  as 

C2{Xk,  (f>k,  fk,  Ik)  =  E  E  { Mn,  n)!"*  -h  [xk{n,  n)  -  Re  [yl{n,  n)e^(2-An+2-Afi+^fc)]]  j 

n=0  n=0  ^  ^ 

Nk-l  Nk-l 

-E  E  {Re2[y*(n,n)e^(2-An+2irAn+^,)]|,  (5.14) 

n=0  n=0 
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where  Re(-)  denotes  the  real  part  and  (•)*  denotes  the  complex  conjugate.  The  estimate  of 
each  element  of  Xk  is  given  by 

x.{n,n)  = 

n=:0,l,---,iVfc-l,  n  =  0,l,---,^ik-l.  (5.15) 

^  ^  £l  — 

Inserting  (5.15)  into  (5.14),  we  obtain  the  NLS  estimates  {(j)k,  fk,  fk}  of  {4)k,  fk,  fk}  by 

equivalently  maximizing  the  following  cost  function: 

Nk-i 

n=0  n=0 

1  ATjt-lJVjt-l 

^  n=0  n=0 

The  estimate  of  (j)k  is  given  by 

i,  =  iangle  Ie'  E'  [i/K".  |  .  (5.17) 

I  >  fk=fkJk=fk 

Inserting  (5.17)  into  (5.16)  and  dropping  out  a  term  that  is  not  related  to  the  frequency 
estimation,  we  finally  obtain  the  estimate  {fk,  fk)  of  [fk,  h}  by  the  following  maximization: 

Nk-lNk-l 

{fk,  /fc}  =  argm^  ^  yl{n,n)e  J2^(2/fcn+2/fcn)  ^  ^5 

fkJk  firrO 

which  can  be  done  by  simply  applying  2-D  FFT  to  y^n,  n)  with  2/*  and  2fk  as  the  frequency 
variables. 

Note  that  Xk{n,n)  usually  does  not  contain  obvious  gaps  due  to  chipping  or  windowing 
in  the  image  domain.  After  the  estimate  Xk{n,n)  of  Xk{n,n)  is  obtained,  the  segments 
of  Xfc(n,  n)  corresponding  to  the  gaps  of  Zk{n,n)  should  be  set  to  zero.  Let  Gk  contain  the 
indices  corresponding  to  where  the  gaps  of  Xk{n,  n)  should  be.  Then  the  relationship  between 
the  indices  in  Gk  and  G  is  linear  and  can  be  computed  with: 

[^0(1)1  ,  (5.19) 
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or 

Gt(i)=  ^G(»)  ,  (5-20) 

where  [•]  and  [-J  denote  the  nearest  integers  towards  positive  infinity  and  minus  infinity, 
respectively,  and  Gk{i)  and  G{i)  denote  the  ith  indices  of  Gk  and  G,  respectively.  We 
compute  the  begining  of  the  gaps  with  (5.19)  and  the  end  of  the  gaps  with  (5.20). 

The  steps  of  the  parameter  estimation  for  a  single  scatterer  can  be  summarized  as  follows: 
Step  1:  Estimate  the  frequency  pair  {fk,  fk}  with  (5.18)  by  using  2-D  FFT  with  zero- 
padding  to  obtain  an  initial  estimate,  which  is  then  refined  by  using  the  FMINS  function  of 
MATLAB. 

Step  2:  Calculate  ^k  according  to  (5.17)  with  {fk,  fk}  replaced  by  [fk,  /*}  obtained  in  Step 

1. 

Step  3:  Determine  Xfc(n,  n),  n  =  0, 1,  -  •  • ,  -  1,  n  =  0, 1,  •  •  • ,  iVfc  -  1,  according  to  (5.15) 

with  {(pk,  fk,  fk}  replaced  by  {^k,  fk,  h}  obtained  in  the  previous  two  steps. 

Step  4:  Compute  the  cross-range  indices  to  form  Gk  and  set  those  columns  of  Xk{n,n) 
corresponding  to  the  indices  in  Gk  to  zero. 

Sine  Function  Based  Real- Valued  Amplitude  Estimation 

When  the  estimate  of  the  2-D  real- valued  amplitude  sequence  Xk{n,n)  for  the  kth  scat¬ 
terer  is  obtained,  we  average  it  over  the  range  dimension  to  obtain 

Vkin)  = ^  Xk{n,n),  n  =  0, 1,  •  •  • ,  JV*,  -  1,  n^Gk-  (5.21) 

n=0 

We  then  approximately  model  the  so-obtained  1-D  real- valued  sequence  with  a  sine  function. 
The  parameter  estimation  of  the  sine  function  is  given  next. 

Let 

yjk(n)  =  afcpjfc(n)  +  efc(n),  n  =  0, 1,  •  •  • ,  Ajt  -  1,  n^Gk-  (5.22) 

where  gk{n)  is  defined  as 

gk{n)  =  sine  [6jt7r(n  -  rt)] .  (5.23) 
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Let  Yfc  and  gk  denote  the  vectors  whose  nth  elements  are  yfc(n)  and  gk{n),  n  ^  Gk, 
respectively.  Then  the  NLS  estimates  of  {ak,  bk,  Tk]  can  be  determined  by  minimizing  the 
following  cost  function 

Criok,  bk,Tk)  =  l|yfc  -  Qkgklf  •  (5-24) 


Minimizing  Cj  in  (5.24)  with  respect  to  ak  gives  the  estimate  ak  of  oik- 


ctk 


EkYk 


iig^ir 


Inserting  (5.25)  into  (5.24),  Cj  in  (5.24)  can  be  simplified  to 

{ylskf 


C^{k,n)  =  llnir  - 

which  can  be  minimized  by  maximizing  its  last  term 


2  ’ 


C9{bk,Tk) 


(yfcgfc)' 

llgfcll" 


(5.25) 


(5.26) 


(5.27) 


The  maximization  of  (5.27)  requires  a  2-D  search  over  the  parameter  space.  We  use  an 
alternating  maximization  procedure  by  iteratively  updating  bk  and  Tk  while  fixing  the  other 
parameter. 

The  initial  estimate  fjt  of  Tk  is  obtained  by  finding  the  peak  position  of  HyitH.  When  we 
search  for  h  with  the  FMIN  function  of  MATLAB,  we  also  need  to  obtain  the  initial  estimate 
bk  of  bk.  We  calculate  the  1-D  FFT  y*/  of  the  data  sequence  y*  and  search  two  frequency 
positions  fi  and  fr  nearest  to  0,  where  <  0  and  fr  >  0,  such  that  |yit/(/i)|  <  2  ly*/(^)l 
and  |yjfc/(/r)|  <  |  |yfc/(0)l-  The  initial  estimate  6*;  of  bk  is  given  by 


^ 

bk  =  fr-  fl- 


(5.28) 


The  algorithm  for  the  NLS  estimation  of  the  parameters  of  the  sine  function  model  in 
(5.22)  is  described  as  follows. 

Step  (1):  Obtain  the  initial  estimates  %  of  Tk  by  using  the  above  mentioned  initialization. 
Step  (2):  Update  bk-  Replace  r*  in  (5.27)  by  fk  obtained  in  the  previous  step  and  search 
for  bk  which  maximizes  Cg. 
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step  (3):  Update  n.  Replace  bk  in  (5.27)  with  bk  obtained  in  Step  (2)  and  search  for  fk 
which  maximizes  Cg. 

Step  (4):  Repeat  Steps  2  and  3  until  “practical  convergence”  which  is  determined  by 
checking  the  relative  change  e  of  the  cost  function  Cg  in  (5.27)  between  two  consecutive 
iterations.  We  use  e  =  10"^  in  our  numerical  examples. 

Step  (5):  Calculate  dfc  using  (5.25). 

5.3.3  Image  Reconstruction 

After  the  estimates  joit,  h,  of  sine  function  parameters  are  obtained,  we  simu¬ 
late  {{^fc(«)}fi=i}^^i  by  ^sing  {djfe,  bk,  according  to  (5.22),  and  repeat{{^fc(n)}fi*i}^^^ 

in  range  to  obtain  the  2-D  real- valued  amplitude  {Xfc}^i,  that  is 

x(n,  n)fc  =  y*.(n),  n  =  0,  !,•••,  74  —  1,  n  =  0, 1,  •  •  •  ,iVfc  —  1,  A:  =  1,2,  •  •  •  ,i^.  (5.29) 

A  A 

With  {Xk}k=i  and  {^/t,  fk,  fk}k^v  fbe  phase  history  data  matrix  of  each  scatterer  is 
simualted  without  gaps.  Then  we  perform  the  2-D  FFT  on  the  simulated  phase  history  data 
matrix  to  obtain  the  enhanced  image  chip  for  each  scatterer.  Each  enhanced  image  chip  is 
put  back  to  its  original  location  in  the  FFT  image  V  defined  in  Section  5.3.1  to  reconstructed 
the  approximate  SAR  image  of  the  target.  If  the  image  chips  of  different  scatterers  overlap, 
we  use  whichever  is  stronger  for  the  overlapping  part. 

5.4  NumericzJ  Results 

We  illustrate  the  SAR  target  feature  extraction  and  imaging  performance  of  QUALE  with 
a  set  of  high  resolution  phase  history  data  generated  by  XPATCH  [19],  a  high  frequency 
electromagnetic  scattering  prediction  code  for  complex  3-D  objects.  A  photo  of  the  slicy 
object  taken  at  45°  azimuth  angle  and  the  corresponding  FFT  SAR  image  from  the  XPATCH 
data  obtained  at  0°  are,  respectively,  shown  in  Figures  5.2(a)  and  (b)  where  a  good  agreement 
between  the  slicy  object  and  the  image  can  be  seen  clearly.  The  XPATCH  data  matrix  we 
use  has  a  resolution  of  0.038  meters  in  both  range  and  cross-range  and  a  size  of  N  =  N  = 


84 


288.  To  demonstrate  the  advantages  of  the  angle  diversity  data  fusion,  we  assume  that  we 
have  three  segments  of  the  XPATCH  phase  history  data  available  to  us  with  each  segment 
consisting  of  continuous  data  columns  from  the  beginning,  middle,  and  end  of  the  original 
data  matrix.  Hence  the  data  are  collected  at  three  diverse  azimuth  angle  locations. 

Now  we  compare  QUALE  with  the  FFT  approaches  as  well  as  CLEAN,  which  has  been 
proved  to  be  very  useful  in  several  applications  including  astronomical  data  analysis,  mi¬ 
crowave  imaging,  spectral  estimation,  and  target  feature  extraction.  When  using  CLEAN, 
we  model  each  scatterer  of  the  target  of  interest  as  a  complex  sinusoid  with  constant  ampli¬ 
tude  and  phase  in  both  range  and  cross-range.  CLEAN  first  estimates  the  parameters  of  the 
strongest  scatterer,  subtracts  it  from  the  original  signal,  and  then  repeats  this  process  for 
the  next  strongest  scatterer  until  the  predetermined  number  of  scatterers  is  reached  or  the 
estimated  amplitude  of  the  current  scatterer  is  small  enough.  Once  the  scatterer  parameters 
are  estimated,  they  are  used  to  simulate  the  data  in  the  missing  gaps.  FFT  is  then  used  to 
form  the  CLEAN  images.  Since  CLEAN  uses  a  point-like  scatterer  model,  a  larger  number 
of  scatterers  have  to  be  assumed  to  model  a  dihedral  corner  reflector.  We  set  a  total  number 
of  scatterers  to  be  30  for  CLEAN. 

The  thresholds  Tr  and  Tcr  used  in  the  isolation  process  of  QUALE  is  2%  and  5%  of 
the  peak  value  in  range  and  cross-range,  respectively.  We  implement  the  scatterer  isolation 
process  of  QUALE  based  on  the  2-D  FFT  of  the  middle  segment  of  the  gapped  data  with 
zero-padding.  In  a  cross-range,  if  the  second  peak  is  greater  than  40%  of  the  largest  one  in 
one  window,  they  will  be  split  into  two  windows.  We  assume  that  the  number  of  scatterers 
is  AT  =  7  for  QUALE. 

We  first  consider  a  case  where  the  gapped  data  set  due  to  angle  diversity  contains  50%  of 
the  original  XPATCH  data  set.  The  original  data  matrix  in  columns  49  through  120  and  169 
through  240  are  set  to  zero  to  form  the  gapped  data  set.  Hence  the  angle  diversity  data  we 
use  in  this  example  has  three  non-zero  data  segments  (each  with  48  columns)  intermitted  by 
two  gaps  (each  with  72  columns).  Figure  5.3(a)  shows  the  windows  used  by  QUALE  to  chip 
out  the  7  scatterers.  Figures  3(b)  and  (c)  show  the  FFT  SAR  images  obtained  by  using  the 
middle  segment  only  and  all  of  the  segments,  respectively.  Compared  with  Figure  5.2(b), 
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Figures  3(b)  and  (c)  exhibit  distinct  energy  leakage  in  cross-range  from  the  mainlobes  of  the 
scatterers,  with  peak  splitting  occurring  in  Figure  5.3(c)  due  to  the  gapped  data.  Figures 
3(d)  shows  the  SAR  image  obtained  by  using  CLEAN.  The  SAR  image  obtained  via  QUALE 
is  presented  in  Figure  5.3(e).  Note  that,  CLEAN  generates  images  with  isolated  points.  The 
line-shaped  features  in  cross-range  caused  by  dihedral-like  scatterers  (see  Figure  5.2)  are 
represented  with  several  point  scatterers.  For  the  trihedral  type  of  corner  reflectors,  CLEAN 
works  very  well  since  the  data  model  is  quite  accurate.  From  Figure  ,5.3 (d)  we  note  that 
QUALE  can  generate  better  images  than  the  FFT  approaches  and  CLEAN  for  dihedral  type 
of  corner  reflectors.  For  this  example,  our  simulations  show  that  the  ratios  between  the 
MATLAB  flops  needed  by  CLEAN  and  QUALE  over  those  needed  by  the  FFT  approaches 
are,  178.73  and  11.59,  respectively.  CLEAN  requires  much  more  computations  than  QUALE 
because  the  former  performs  2-D  FFT  on  the  entire  data  matrix  while  QUALE  works  with 
the  scatterer  image  chips  with  much  smaller  dimensions.  Also,  the  NLS  approach  used  to 
estimate  the  sine  function  parameters  in  QUALE  causes  little  additional  burden  to  the  total 
amount  of  computations  since  the  approach  convergences  quickly  within  a  few  iterations. 

Consider  next  an  example  with  much  larger  gaps.  We  set  the  original  XPATCH  data 
matrix  in  columns  11  through  139  and  150  through  278  to  zero  to  form  the  angle  data 
diversity  data  set.  Three  data  segments  each  with  10  columns  are  separated  by  two  gaps 
each  with  129  columns.  That  is,  the  total  angle  diversity  data  set  contains  only  10%  of 
the  original  data  set.  The  cross-range  resolution  of  each  segment  is  288  x  0.038/10  «  1.094 
meters.  The  windows  used  by  QUALE  to  chip  out  the  7  scatterers  is  shown  in  Figure  5.4(a). 
The  SAR  images  obtained  by  using  the  FFT  approaches,  CLEAN,  and  QUALE  are  shown 
in  Figures  5.4(b)  through  (e),  which  correspond  to  the  images  in  Figures  5.3(b)  through  (e), 
respectively.  A  comparison  between  Figures  5.3  and  5.4  shows  that  the  FFT  SAR  images 
become  worse  due  to  the  larger  data  gaps  with  more  mainlobe  energy  leaked  in  cross-range 
and  more  severe  peak  splitting.  While  both  FFT  approaches  and  CLEAN  need  the  same 
amounts  of  computations  as  in  the  previous  example,  QUALE  requires  approximately  40% 
more  computations  than  in  the  previous  example  due  to  larger  image  chips.  Hence  QUALE 
is  still  computationally  much  more  efficient  than  CLEAN. 


86 


5.5  Conclusions 


We  have  presented  a  quasi-parametric  estimation  algorithm,  referred  to  as  QUALE,  for 
approximate  SAR  target  feature  extraction  and  imaging  with  angle  diversity  based  on  a 
flexible  data  model,  which  describes  each  target  scatterer  as  a  2-D  complex  sequence  with 
arbitrary  amplitude  and  constant  phase  in  range  and  cross-range.  This  data  model  is  es¬ 
sentially  quasi-parametric  because  of  the  arbitrary  amplitude  assumption.  QUALE  first 
estimates  the  model  parameters  that  include,  for  each  scatterer,  a  2-D  arbitrary  real- valued 
amplitude  sequence,  a  constant  phase,  and  the  locations  in  range  and  cross-range.  It  then 
averages  the  estimated  2-D  real-valued  amplitude  sequence  over  the  range  dimension,  ap¬ 
proximately  models  the  so-obtained  1-D  sequence  with  a  sine  function,  and  estimates  the 
sine  function  parameters  by  minimizing  a  NLS  cost  function.  Finally,  the  2-D  SAR  image 
is  reconstructed  by  using  the  estimated  features.  We  have  shown  with  XPATCH  simulated 
examples  that  QUAI.E  outperforms  the  FFT  approaches  and  the  well-known  CLEAN  algo¬ 
rithm  for  a  target  containing  both  dihedrals  and  trihedrals  even  in  the  presence  of  very  large 
data  gaps.  QUALE  is  also  computationally  much  more  efficient  than  CLEAN. 
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(a) 


(b) 

Figure  5.1;  1-D  FFTs  of  continuous  data  and  gapped  data  due  to  angle  diversity.  The 
continuous  data  used  is  an  all-one  sequence  with  length  32,  and  the  gapped  data  due  to 
angle  diversity  is  formed  by  setting  the  middle  50%  of  the  continuous  data  to  zero,  (a)  FFT 
of  continuous  data,  (b)  FFT  of  gapped  data  due  to  angle  diversity. 
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(a) 


Figure  5.2:  (a)  Target  photo  taken  at  45°  azimuth  angle,  (b)  2-D  FFT  SAR  image  obtained 
at  0°  azimuth  angle  using  the  entire  XPATCH  data  with  a  resolution  0.038  meters  in  both 
range  and  cross-range. 
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(e) 

Figure  5.3:  Comparison  of  SAR  images  obtained  via  different  algorithms  for  the  XPATCH 
data  with  two  gaps  in  cross-range  (only  50%  of  the  data  are  used  with  cross-range  resolution 
of  each  segment  being  288  x  0.038/48  w  0.228  meters),  (a)  The  windows  used  by  QUALE  to 
chip  out  the  7  scatterers,  (b)  2-D  FFT  SAR  image  from  the  middle  segment  of  the  gapped 
XPATCH  data  with  zero-padding,  (c)  2-D  FFT  SAR  image  from  the  entire  gapped  XPATCH 
data,  (d)  CLEAN  SAR  image,  and  (d)  QUALE  SAR  image. 
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(d)  (e) 

Figure  5.4:  Comparison  of  SAR  images  obtained  via  different  algorithms  for  the  XPATCH 
data  with  two  gaps  in  cross- range  (only  10%  of  the  data  are  used  with  cross-range  resolution 
of  each  segment  being  288  x  0.038/10  «  1.094  meters),  (a)  The  windows  used  by  QUALE  to 
chip  out  the  7  scatterers,  (b)  2-D  FFT  SAR  image  from  the  middle  segment  of  the  gapped 
XPATCH  data  with  zero-padding,  (c)  2-D  FFT  SAR  image  from  the  entire  gapped  XPATCH 
data,  (d)  CLEAN  SAR  image,  and  (d)  QUALE  SAR  image. 
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